FV3 Bundle
nwp_nudge_nlm.F90
Go to the documentation of this file.
1 module nwp_nudge_nlm_mod
2 
3  use fv_arrays_nlm_mod, only: real4, real8, fvprc
4  use fv_control_nlm_mod, only: npx, npy
5  use fv_grid_utils_nlm_mod, only: i_sst, j_sst, sst_ncep
6  use fv_grid_utils_nlm_mod, only: vlon, vlat, sina_u, sina_v, da_min, great_circle_dist, ks, intp_great_circle
7  use fv_grid_tools_nlm_mod, only: agrid, dx, dy, rdxc, rdyc, rarea, area
8  use fv_diagnostics_nlm_mod,only: prt_maxmin, fv_time
9  use tp_core_nlm_mod, only: copy_corners
10  use fv_mapz_nlm_mod, only: mappm
11  use fv_mp_nlm_mod, only: is,js,ie,je, isd,jsd,ied,jed, gid, masterproc, domain, mp_reduce_sum
12  use fv_timing_nlm_mod, only: timing_on, timing_off
13  use constants_mod, only: pi, grav, rdgas, cp_air, kappa, radius
14  use time_manager_mod, only: time_type, get_time, get_date
15  use mpp_mod, only: mpp_error, fatal, stdlog
16  use fms_mod, only: write_version_number, open_namelist_file, &
17  check_nml_error, file_exist, close_file, &
18  read_data, field_exist
19  use fms_io_mod, only: field_size
20  use mpp_domains_mod, only: mpp_update_domains
21 
22  implicit none
23  private
24 
25  character(len=128) :: version = ''
26  character(len=128) :: tagname = ''
27 
28  public do_nwp_nudge, nwp_nudge_init, nwp_nudge_end, breed_slp_inline
29 
30  integer im ! Data x-dimension
31  integer jm ! Data y-dimension
32  integer km ! Data z-dimension
33  real(FVPRC), allocatable:: ak0(:), bk0(:)
34  real(FVPRC), allocatable:: lat(:), lon(:)
35 
36  logical :: module_is_initialized = .false.
37  logical :: master
38  logical :: no_obs
39  real(FVPRC) :: deg2rad, rad2deg
40  real(FVPRC) :: time_nudge = 0.
41  integer :: time_interval = 6*3600 ! dataset time interval (seconds)
42  integer, parameter :: nfile_max = 125
43  integer :: nfile = 1
44 
45  integer :: k_breed = 0
46  integer :: k_trop = 0
47  real(FVPRC) :: p_trop = 300.e2
48 
49  real(FVPRC), allocatable:: s2c(:,:,:)
50  integer, allocatable:: id1(:,:), id2(:,:), jdc(:,:)
51  real(FVPRC), allocatable :: ps_dat(:,:,:)
52  real*4, allocatable, dimension(:,:,:,:):: u_dat, v_dat, t_dat, q_dat
53  real(FVPRC), allocatable:: gz0(:,:)
54 
55 ! Namelist variables:
56  character(len=128):: file_names(nfile_max)
57  character(len=128):: track_file_name
58  integer :: nfile_total = 0 ! =5 for 1-day (if datasets are 6-hr apart)
59  real(FVPRC) :: p_wvp = 100.e2 ! cutoff level for specific humidity nudging
60  integer :: kord_data = 8
61 
62  logical :: tc_mask = .false.
63  logical :: strong_mask = .true.
64  logical :: ibtrack = .false.
65  logical :: nudge_debug = .false.
66  logical :: nudge_t = .false.
67  logical :: nudge_q = .false.
68  logical :: nudge_winds = .true.
69  logical :: nudge_virt = .false.
70  logical :: nudge_hght = .false.
71  logical :: nudge_tpw = .false. ! nudge total precipitable water
72  logical :: time_varying = .true.
73  logical :: time_track = .false.
74 
75 
76 ! Nudging time-scales (seconds): note, however, the effective time-scale is 2X smaller (stronger) due
77 ! to the use of the time-varying weighting factor
78  real(FVPRC) :: tau_q = 86400. ! 1-day
79  real(FVPRC) :: tau_tpw = 86400. ! 1-day
80  real(FVPRC) :: tau_winds = 21600. ! 6-hr
81  real(FVPRC) :: tau_t = 86400.
82  real(FVPRC) :: tau_virt = 86400.
83  real(FVPRC) :: tau_hght = 86400.
84 
85  real(FVPRC) :: q_min = 1.e-8
86  real(FVPRC) :: q_rat
87 
88  integer :: nf_uv = 0
89  integer :: nf_t = 2
90 
91 ! starting layer (top layer is sponge layer and is skipped)
92  integer :: kstart = 2
93 
94 ! skip "kbot" layers
95  integer :: kbot_winds = 0
96  integer :: kbot_t = 0
97  integer :: kbot_q = 1
98 
99 !-- Tropical cyclones --------------------------------------------------------------------
100 
101 ! track dataset: 'INPUT/tropical_cyclones.txt'
102 
103  logical :: breed_vortex = .false.
104  real(FVPRC) :: tau_vortex = 600.
105 
106  real(FVPRC) :: slp_env = 101010. ! storm environment pressure (pa)
107  real(FVPRC) :: r_min = 200.e3
108  real(FVPRC) :: r_inc = 25.e3
109  real(FVPRC), parameter:: del_r = 50.e3
110  integer:: year_track_data
111  integer, parameter:: max_storm = 140 ! max # of storms to process
112  integer, parameter:: nobs_max = 125 ! Max # of observations per storm
113 
114  integer :: nstorms = 0
115  integer :: nobs_tc(max_storm)
116  real*4:: x_obs(nobs_max,max_storm) ! longitude in degrees
117  real*4:: y_obs(nobs_max,max_storm) ! latitude in degrees
118  real*4:: mslp_obs(nobs_max,max_storm) ! observed SLP in mb
119  real*4:: mslp_out(nobs_max,max_storm) ! outer ring SLP in mb
120  real*4:: rad_out(nobs_max,max_storm) ! outer ring radius in meters
121  real*4:: time_tc(nobs_max,max_storm) ! start time of the track
122 !------------------------------------------------------------------------------------------
123 
124  namelist /nwp_nudge_nml/ nudge_virt, nudge_hght, nudge_t, nudge_q, nudge_winds, nudge_tpw, &
125  tau_winds, tau_t, tau_q, tau_virt, tau_hght, kstart, kbot_winds, &
126  k_breed, k_trop, p_trop, kord_data, tc_mask, nudge_debug, nf_t, &
127  nf_uv, breed_vortex, tau_vortex, tau_tpw, strong_mask, &
128  kbot_t, kbot_q, p_wvp, time_varying, time_track, time_interval, &
129  r_min, r_inc, ibtrack, track_file_name, file_names
130 
131  contains
132 
133 
134  subroutine do_nwp_nudge ( Time, dt, npz, ps_dt, u_dt, v_dt, t_dt, q_dt, zvir, &
135  ak, bk, ts, ps, delp, ua, va, pt, nwat, q, phis )
136 
137  type(time_type), intent(in):: time
138  integer, intent(in):: npz ! vertical dimension
139  integer, intent(in):: nwat
140  real(FVPRC), intent(in):: dt
141  real(FVPRC), intent(in):: zvir
142  real(FVPRC), intent(in ), dimension(npz+1):: ak, bk
143  real(FVPRC), intent(in ), dimension(isd:ied,jsd:jed ):: phis
144  real(FVPRC), intent(inout), dimension(isd:ied,jsd:jed,npz):: pt, ua, va, delp
145  real(FVPRC), intent(inout):: q(isd:ied,jsd:jed,npz,nwat)
146  real(FVPRC), intent(inout), dimension(isd:ied,jsd:jed):: ps
147 ! Accumulated tendencies
148  real(FVPRC), intent(inout), dimension(isd:ied,jsd:jed,npz):: u_dt, v_dt
149  real(FVPRC), intent(out):: t_dt(is:ie,js:je,npz)
150  real(FVPRC), intent(out):: q_dt(is:ie,js:je,npz)
151  real(FVPRC), intent(out), dimension(is:ie,js:je):: ps_dt, ts
152 ! local:
153  real(FVPRC):: tpw_dat(is:ie,js:je)
154  real(FVPRC):: tpw_mod(is:ie,js:je)
155  real(FVPRC):: mask(is:ie,js:je)
156  real(FVPRC):: gz_int(is:ie,js:je), gz(is:ie,npz+1), peln(is:ie,npz+1), pk(is:ie,npz+1)
157  real(FVPRC):: pe1(is:ie)
158  real(FVPRC):: pkz, ptmp
159  real(FVPRC), allocatable :: ps_obs(:,:)
160  real(FVPRC), allocatable, dimension(:,:,:):: u_obs, v_obs, t_obs, q_obs
161  integer :: seconds, days
162  integer :: i,j,k, iq, kht
163  real(FVPRC) :: factor
164  real(FVPRC) :: dbk, rdt, press(npz), profile(npz), prof_t(npz), prof_q(npz), du, dv
165 
166 
167  if ( .not. module_is_initialized ) then
168  call mpp_error(fatal,'==> Error from do_nwp_nudge: module not initialized')
169  endif
170 
171  call get_time (time, seconds, days)
172 
173  do j=js,je
174  do i=is,ie
175  mask(i,j) = 1.
176  enddo
177  enddo
178  if ( tc_mask ) call get_tc_mask(time, mask)
179 
180 ! The following profile is suitable only for NWP purposes; if the analysis has a good representation
181 ! of the strat-meso-sphere the profile for upper layers should be changed.
182 
183  profile(:) = 1.
184  do k=1,npz
185  press(k) = 0.5*(ak(k) + ak(k+1)) + 0.5*(bk(k)+bk(k+1))*1.e5
186  if ( press(k) < 30.e2 ) then
187  profile(k) = max(0.01, press(k)/30.e2)
188  endif
189  enddo
190  profile(1) = 0.
191 
192 ! Thermodynamics:
193  prof_t(:) = 1.
194  do k=1,npz
195  if ( press(k) < 30.e2 ) then
196  prof_t(k) = max(0.01, press(k)/30.e2)
197  endif
198  enddo
199  prof_t(1) = 0.
200 
201 ! Water vapor:
202  prof_q(:) = 1.
203  do k=1,npz
204  if ( press(k) < 300.e2 ) then
205  prof_q(k) = max(0., press(k)/300.e2)
206  endif
207  enddo
208  prof_q(1) = 0.
209 
210 ! Height
211  if ( k_trop == 0 ) then
212  k_trop = 2
213  do k=2,npz-1
214  ptmp = ak(k+1) + bk(k+1)*1.e5
215  if ( ptmp > p_trop ) then
216  k_trop = k
217  exit
218  endif
219  enddo
220  endif
221 
222  if ( time_varying ) then
223  factor = 1. + cos(real(mod(seconds,time_interval))/real(time_interval)*2.*pi)
224  factor = max(1.e-5, factor)
225  else
226  factor = 1.
227  endif
228 
229  allocate (ps_obs(is:ie,js:je) )
230  allocate ( t_obs(is:ie,js:je,npz) )
231  allocate ( q_obs(is:ie,js:je,npz) )
232 
233  if ( nudge_winds ) then
234  allocate (u_obs(is:ie,js:je,npz) )
235  allocate (v_obs(is:ie,js:je,npz) )
236  endif
237 
238 
239  call get_obs(time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, u_obs, v_obs, t_obs, q_obs, &
240  tpw_dat, phis, gz_int, npz)
241 
242  if ( no_obs ) return
243 
244  ps_dt(:,:) = 0.
245 
246  if ( nudge_winds ) then
247 
248 ! Compute tendencies:
249  rdt = 1. / (tau_winds/factor + dt)
250  do k=kstart, npz - kbot_winds
251  do j=js,je
252  do i=is,ie
253  u_obs(i,j,k) = profile(k)*(u_obs(i,j,k)-ua(i,j,k))*rdt
254  v_obs(i,j,k) = profile(k)*(v_obs(i,j,k)-va(i,j,k))*rdt
255  enddo
256  enddo
257  enddo
258 
259  if ( nf_uv>0 ) call del2_uv(u_obs, v_obs, 0.20, npz, nf_uv)
260 
261  do k=kstart, npz - kbot_winds
262  do j=js,je
263  do i=is,ie
264 ! Apply TC mask
265  u_obs(i,j,k) = u_obs(i,j,k) * mask(i,j)
266  v_obs(i,j,k) = v_obs(i,j,k) * mask(i,j)
267 !
268  u_dt(i,j,k) = u_dt(i,j,k) + u_obs(i,j,k)
269  v_dt(i,j,k) = v_dt(i,j,k) + v_obs(i,j,k)
270  ua(i,j,k) = ua(i,j,k) + u_obs(i,j,k)*dt
271  va(i,j,k) = va(i,j,k) + v_obs(i,j,k)*dt
272  enddo
273  enddo
274  enddo
275  deallocate ( u_obs )
276  deallocate ( v_obs )
277  endif
278 
279  if ( nudge_t .or. nudge_virt ) then
280  if(nudge_debug) call prt_maxmin('T_obs', t_obs, is, ie, js, je, 0, npz, 1., master)
281  endif
282 
283 !---------------------- temp -----------
284  if ( nudge_virt .and. nudge_hght ) then
285  tau_virt = max(tau_hght, tau_virt)
286  kht = k_trop
287  else
288  kht = npz-kbot_t
289  endif
290 !---------------------- temp -----------
291 
292  t_dt(:,:,:) = 0.
293 
294  if ( nudge_hght ) then
295  if(nudge_debug) call prt_maxmin('H_int', gz_int, is, ie, js, je, 0, 1, 1./grav, master)
296 
297  rdt = dt / (tau_hght/factor + dt)
298 
299  do j=js,je
300 
301  do i=is,ie
302  pe1(i) = ak(1)
303  peln(i,1) = log(pe1(i))
304  pk(i,1) = pe1(i)**kappa
305  enddo
306  do k=2, npz+1
307  do i=is,ie
308  pe1(i) = pe1(i) + delp(i,j,k-1)
309  peln(i,k) = log(pe1(i))
310  pk(i,k) = pe1(i)**kappa
311  enddo
312  enddo
313 
314  do i=is,ie
315  gz(i,npz+1) = phis(i,j)
316  enddo
317  do i=is,ie
318  do k=npz,k_trop+1,-1
319  gz(i,k) = gz(i,k+1) + rdgas*pt(i,j,k)*(1.+zvir*q(i,j,k,1))*(peln(i,k+1)-peln(i,k))
320  enddo
321  enddo
322 
323  do i=is,ie
324  do k=k_trop+1,npz
325 ! Add constant "virtual potential temperature" increment to correct height at p_interface
326  pkz = (pk(i,k+1)-pk(i,k))/(kappa*(peln(i,k+1)-peln(i,k)))
327  pt(i,j,k) = pt(i,j,k) + mask(i,j)*rdt*pkz*(gz_int(i,j)-gz(i,k_trop+1)) / &
328  (cp_air*(1.+zvir*q(i,j,k,1))*(pk(i,npz+1)-pk(i,k_trop+1)))
329  enddo
330  enddo
331  enddo ! j-loop
332  endif
333 
334  if ( nudge_t ) then
335  rdt = 1./(tau_t/factor + dt)
336  do k=kstart, kht
337  do j=js,je
338  do i=is,ie
339  t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)-pt(i,j,k))*rdt
340  enddo
341  enddo
342  enddo
343  elseif ( nudge_virt ) then
344  rdt = 1./(tau_virt/factor + dt)
345  do k=kstart, kht
346  do j=js,je
347  do i=is,ie
348  t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)/(1.+zvir*q(i,j,k,1))-pt(i,j,k))*rdt
349  enddo
350  enddo
351  enddo
352  endif
353 
354  deallocate ( t_obs )
355 
356  if ( nudge_t .or. nudge_virt ) then
357 ! Filter t_dt here:
358  if ( nf_t>0 ) call del2_scalar(t_dt, 0.20, npz, nf_t)
359 
360  do k=kstart, kht
361  do j=js,je
362  do i=is,ie
363  pt(i,j,k) = pt(i,j,k) + t_dt(i,j,k)*dt*mask(i,j)
364  enddo
365  enddo
366  enddo
367  endif
368 
369  q_dt(:,:,:) = 0.
370  if ( nudge_q ) then
371  rdt = 1./(tau_q/factor + dt)
372  do k=kstart, npz - kbot_q
373  if ( press(k) > p_wvp ) then
374  do iq=2,nwat
375  do j=js,je
376  do i=is,ie
377  q(i,j,k,iq) = q(i,j,k,iq)*delp(i,j,k)
378  enddo
379  enddo
380  enddo
381 ! Specific humidity:
382  do j=js,je
383  do i=is,ie
384  delp(i,j,k) = delp(i,j,k)*(1.-q(i,j,k,1))
385  q_dt(i,j,k) = prof_q(k)*(max(q_min,q_obs(i,j,k))-q(i,j,k,1))*rdt*mask(i,j)
386  q(i,j,k,1) = q(i,j,k,1) + q_dt(i,j,k)*dt
387  delp(i,j,k) = delp(i,j,k)/(1.-q(i,j,k,1))
388  enddo
389  enddo
390  do iq=2,nwat
391  do j=js,je
392  do i=is,ie
393  q(i,j,k,iq) = q(i,j,k,iq)/delp(i,j,k)
394  enddo
395  enddo
396  enddo
397  endif
398  enddo
399  elseif ( nudge_tpw ) then
400 ! Compute tpw_model
401  tpw_mod(:,:) = 0.
402  do k=1,npz
403  do j=js,je
404  do i=is,ie
405  tpw_mod(i,j) = tpw_mod(i,j) + q(i,j,k,1)*delp(i,j,k)
406  enddo
407  enddo
408  enddo
409 
410  do j=js,je
411  do i=is,ie
412  tpw_dat = tpw_dat(i,j) / max(tpw_mod(i,j), q_min)
413  enddo
414  enddo
415  if(nudge_debug) call prt_maxmin('TPW_rat', tpw_dat, is, ie, js, je, 0, 1, 1., master)
416 
417  do k=1,npz
418 
419  do iq=2,nwat
420  do j=js,je
421  do i=is,ie
422  q(i,j,k,iq) = q(i,j,k,iq)*delp(i,j,k)
423  enddo
424  enddo
425  enddo
426 
427  do j=js,je
428  do i=is,ie
429  delp(i,j,k) = delp(i,j,k)*(1.-q(i,j,k,1))
430  q_rat = max(0.25, min(tpw_dat(i,j), 4.))
431  q(i,j,k,1) = q(i,j,k,1)*(tau_tpw/factor + mask(i,j)*dt*q_rat)/(tau_tpw/factor + mask(i,j)*dt)
432  delp(i,j,k) = delp(i,j,k)/(1.-q(i,j,k,1))
433  enddo
434  enddo
435 
436  do iq=2,nwat
437  do j=js,je
438  do i=is,ie
439  q(i,j,k,iq) = q(i,j,k,iq)/delp(i,j,k)
440  enddo
441  enddo
442  enddo
443 
444  enddo ! k-loop
445  endif
446 
447  deallocate ( q_obs )
448  deallocate ( ps_obs )
449 
450 
451 
452  if ( breed_vortex ) &
453  call breed_slp(time, dt, npz, ak, bk, ps, phis, delp, ua, va, u_dt, v_dt, pt, q, nwat, zvir)
454 
455  end subroutine do_nwp_nudge
456 
457 
458  subroutine get_obs(Time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, u_obs, v_obs, t_obs, q_obs, &
459  tpw_dat, phis, gz_int, npz)
460  type(time_type), intent(in):: Time
461  integer, intent(in):: npz ! vertical dimension
462  real(FVPRC), intent(in):: zvir
463  real(FVPRC), intent(in):: dt
464  real(FVPRC), intent(in), dimension(npz+1):: ak, bk
465  real(FVPRC), intent(in), dimension(isd:ied,jsd:jed):: phis
466  real(FVPRC), intent(in), dimension(isd:ied,jsd:jed,npz):: delp
467  real(FVPRC), intent(inout), dimension(isd:ied,jsd:jed):: ps
468  real(FVPRC), intent(out), dimension(is:ie,js:je):: ts, ps_obs
469  real(FVPRC), intent(out), dimension(is:ie,js:je,npz):: u_obs, v_obs, t_obs, q_obs
470  real(FVPRC), intent(out):: gz_int(is:ie,js:je)
471  real(FVPRC), intent(out):: tpw_dat(is:ie,js:je)
472 ! local:
473  real*4, allocatable:: ut(:,:,:), vt(:,:,:)
474  real(FVPRC), dimension(is:ie,js:je):: h1, h2
475  integer :: seconds, days
476  integer :: i,j,k
477  real(FVPRC) :: alpha, beta
478 
479  call get_time (time, seconds, days)
480 
481  seconds = seconds - nint(dt)
482 
483 ! Data must be "time_interval" (hr) apart; keep two time levels in memory
484 
485  no_obs = .false.
486 
487  if ( mod(seconds, time_interval) == 0 ) then
488 
489  if ( nfile > nfile_total ) then
490  no_obs = .true.
491  return ! free-running mode
492  else
493 
494  ps_dat(:,:,1) = ps_dat(:,:,2)
495  if ( nudge_winds ) then
496  u_dat(:,:,:,1) = u_dat(:,:,:,2)
497  v_dat(:,:,:,1) = v_dat(:,:,:,2)
498  endif
499  t_dat(:,:,:,1) = t_dat(:,:,:,2)
500  q_dat(:,:,:,1) = q_dat(:,:,:,2)
501 
502 !---------------
503 ! Read next data
504 !---------------
505  call get_ncep_analysis ( ps_dat(:,:,2), u_dat(:,:,:,2), v_dat(:,:,:,2), &
506  t_dat(:,:,:,2), q_dat(:,:,:,2), zvir, &
507  ts, nfile, file_names(nfile) )
508  time_nudge = dt
509  endif
510  else
511  time_nudge = time_nudge + dt
512  endif
513 
514 !--------------------
515 ! Time interpolation:
516 !--------------------
517 
518  beta = time_nudge / real(time_interval)
519 
520  if ( beta < 0. .or. beta > (1.+1.e-7) ) then
521  call mpp_error(fatal,'==> Error from get_obs:data out of range')
522  endif
523 
524  alpha = 1. - beta
525 
526 ! Warning: ps_data is not adjusted for the differences in terrain yet
527  ps_obs(:,:) = alpha*ps_dat(:,:,1) + beta*ps_dat(:,:,2)
528 
529  allocate ( ut(is:ie,js:je,npz) )
530  allocate ( vt(is:ie,js:je,npz) )
531 
532  if ( nudge_winds ) then
533 
534  call remap_uv(npz, ak, bk, ps(is:ie,js:je), delp, ut, vt, &
535  km, ps_dat(is:ie,js:je,1), u_dat(:,:,:,1), v_dat(:,:,:,1) )
536 
537  u_obs(:,:,:) = alpha*ut(:,:,:)
538  v_obs(:,:,:) = alpha*vt(:,:,:)
539 
540  call remap_uv(npz, ak, bk, ps(is:ie,js:je), delp, ut, vt, &
541  km, ps_dat(is:ie,js:je,2), u_dat(:,:,:,2), v_dat(:,:,:,2) )
542 
543  u_obs(:,:,:) = u_obs(:,:,:) + beta*ut(:,:,:)
544  v_obs(:,:,:) = v_obs(:,:,:) + beta*vt(:,:,:)
545  endif
546 
547  if ( nudge_t .or. nudge_virt .or. nudge_q .or. nudge_tpw ) then
548 
549  call remap_tq(npz, ak, bk, ps(is:ie,js:je), delp, ut, vt, &
550  km, ps_dat(is:ie,js:je,1), t_dat(:,:,:,1), q_dat(:,:,:,1), zvir)
551 
552  t_obs(:,:,:) = alpha*ut(:,:,:)
553  q_obs(:,:,:) = alpha*vt(:,:,:)
554 
555  call remap_tq(npz, ak, bk, ps(is:ie,js:je), delp, ut, vt, &
556  km, ps_dat(is:ie,js:je,2), t_dat(:,:,:,2), q_dat(:,:,:,2), zvir)
557 
558  t_obs(:,:,:) = t_obs(:,:,:) + beta*ut(:,:,:)
559  q_obs(:,:,:) = q_obs(:,:,:) + beta*vt(:,:,:)
560 
561  do j=js,je
562  do i=is,ie
563  tpw_dat(i,j) = 0.
564  enddo
565  enddo
566  if ( nudge_tpw ) then
567  do k=1,km
568  do j=js,je
569  do i=is,ie
570  tpw_dat(i,j) = tpw_dat(i,j) + q_obs(i,j,k) * &
571  ( ak0(k+1)-ak0(k) + (bk0(k+1)-bk0(k))*ps_obs(i,j) )
572  enddo
573  enddo
574  enddo
575  endif
576  endif
577 
578  if ( nudge_hght ) then
579  call get_int_hght(h1, npz, ak, bk, ps(is:ie,js:je), delp, ps_dat(is:ie,js:je,1), t_dat(:,:,:,1))
580 ! if(nudge_debug) call prt_maxmin('H_1', h1, is, ie, js, je, 0, 1, 1./grav, master)
581 
582  call get_int_hght(h2, npz, ak, bk, ps(is:ie,js:je), delp, ps_dat(is:ie,js:je,2), t_dat(:,:,:,2))
583 ! if(nudge_debug) call prt_maxmin('H_2', h2, is, ie, js, je, 0, 1, 1./grav, master)
584 
585  gz_int(:,:) = alpha*h1(:,:) + beta*h2(:,:)
586  endif
587 
588  deallocate ( ut )
589  deallocate ( vt )
590 
591  end subroutine get_obs
592 
593 
594  subroutine nwp_nudge_init(npz, zvir, ak, bk, ts, phis)
595  integer, intent(in):: npz ! vertical dimension
596  real(FVPRC), intent(in):: zvir
597  real(FVPRC), intent(in), dimension(isd:ied,jsd:jed):: phis
598  real(FVPRC), intent(in), dimension(npz+1):: ak, bk
599  real(FVPRC), intent(out), dimension(is:ie,js:je):: ts
600  logical found
601  integer tsize(4)
602  integer :: i, j, unit, io, ierr, nt, k
603 
604  master = gid==masterproc
605 
606  deg2rad = pi/180.
607  rad2deg = 180./pi
608 
609  do nt=1,nfile_max
610  file_names(nt) = "No_File_specified"
611  enddo
612 
613  track_file_name = "No_File_specified"
614 
615  if( file_exist( 'input.nml' ) ) then
616  unit = open_namelist_file()
617  io = 1
618  do while ( io .ne. 0 )
619  read( unit, nml = nwp_nudge_nml, iostat = io, end = 10 )
620  ierr = check_nml_error(io,'nwp_nudge_nml')
621  end do
622 10 call close_file ( unit )
623  end if
624  call write_version_number (version, tagname)
625  if ( master ) then
626  write( stdlog(), nml = nwp_nudge_nml )
627  write(*,*) 'NWP nudging initialized.'
628  endif
629 
630  if ( nudge_virt ) then
631  nudge_t = .false.
632 ! nudge_q = .false.
633  endif
634 
635  if ( nudge_t ) nudge_virt = .false.
636 
637  do nt=1,nfile_max
638  if ( file_names(nt) == "No_File_specified" ) then
639  nfile_total = nt - 1
640  if(master) write(*,*) 'Total of NCEP files specified=', nfile_total
641  exit
642  endif
643  enddo
644 
645 
646 ! Initialize remapping coefficients:
647 
648  call field_size(file_names(1), 'T', tsize, field_found=found)
649 
650  if ( found ) then
651  im = tsize(1); jm = tsize(2); km = tsize(3)
652  if(master) write(*,*) 'NCEP analysis dimensions:', tsize
653  else
654  call mpp_error(fatal,'==> Error from get_ncep_analysis: T field not found')
655  endif
656 
657  allocate ( s2c(is:ie,js:je,4) )
658  allocate ( id1(is:ie,js:je) )
659  allocate ( id2(is:ie,js:je) )
660  allocate ( jdc(is:ie,js:je) )
661 
662  allocate ( lon(im) )
663  allocate ( lat(jm) )
664 
665  call read_data (file_names(1), 'LAT', lat, no_domain=.true.)
666  call read_data (file_names(1), 'LON', lon, no_domain=.true.)
667 
668 ! Convert to radian
669  do i=1,im
670  lon(i) = lon(i) * deg2rad ! lon(1) = 0.
671  enddo
672  do j=1,jm
673  lat(j) = lat(j) * deg2rad
674  enddo
675 
676  allocate ( ak0(km+1) )
677  allocate ( bk0(km+1) )
678 
679  call read_data (file_names(1), 'hyai', ak0, no_domain=.true.)
680  call read_data (file_names(1), 'hybi', bk0, no_domain=.true.)
681 
682 ! Note: definition of NCEP hybrid is p(k) = a(k)*1.E5 + b(k)*ps
683  ak0(:) = ak0(:) * 1.e5
684 
685 ! Limiter to prevent NAN at top during remapping
686  ak0(1) = max(1.e-8, ak0(1))
687 
688  if ( master ) then
689  do k=1,npz
690  write(*,*) k, 0.5*(ak(k)+ak(k+1))+0.5*(bk(k)+bk(k+1))*1.e5, 'del-B=', bk(k+1)-bk(k)
691  enddo
692  endif
693 
694  if ( k_breed==0 ) k_breed = ks
695 ! k_breed = ks
696 
697  call slp_obs_init
698 
699 !-----------------------------------------------------------
700 ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
701 !-----------------------------------------------------------
702  call remap_coef
703 
704  allocate ( gz0(is:ie,js:je) )
705  allocate (ps_dat(is:ie,js:je,2) )
706  allocate ( u_dat(is:ie,js:je,km,2) )
707  allocate ( v_dat(is:ie,js:je,km,2) )
708  allocate ( t_dat(is:ie,js:je,km,2) )
709  allocate ( q_dat(is:ie,js:je,km,2) )
710 
711 
712 ! Get first dataset
713  nt = 2
714  call get_ncep_analysis ( ps_dat(:,:,nt), u_dat(:,:,:,nt), v_dat(:,:,:,nt), &
715  t_dat(:,:,:,nt), q_dat(:,:,:,nt), zvir, &
716  ts, nfile, file_names(nfile) )
717 
718 
719  module_is_initialized = .true.
720 
721  end subroutine nwp_nudge_init
722 
723 
724  subroutine get_ncep_analysis ( ps, u, v, t, q, zvir, ts, nfile, fname )
725  real(FVPRC), intent(in):: zvir
726  character(len=128), intent(in):: fname
727  integer, intent(inout):: nfile
728 !
729  real(FVPRC), intent(out), dimension(is:ie,js:je):: ts
730  real(FVPRC), intent(out), dimension(is:ie,js:je):: ps
731  real*4, intent(out), dimension(is:ie,js:je,km):: u, v, t, q
732 ! local:
733  real(FVPRC), allocatable:: oro(:,:), wk2(:,:), wk3(:,:,:)
734  real(FVPRC) tmean
735  integer:: i, j, k, npt
736  integer:: i1, i2, j1
737  logical found
738  logical:: read_ts = .true.
739  logical:: land_ts = .false.
740 
741  if( .not. file_exist(fname) ) then
742  call mpp_error(fatal,'==> Error from get_ncep_analysis: file not found')
743  else
744  if(master) write(*,*) 'Reading NCEP anlysis file:', fname
745  endif
746 
747 !----------------------------------
748 ! remap surface pressure and height:
749 !----------------------------------
750  allocate ( wk2(im,jm) )
751  call read_data (fname, 'PS', wk2, no_domain=.true.)
752  if(gid==0) call pmaxmin( 'PS_ncep', wk2, im, jm, 0.01)
753 
754  do j=js,je
755  do i=is,ie
756  i1 = id1(i,j)
757  i2 = id2(i,j)
758  j1 = jdc(i,j)
759  ps(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
760  s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
761  enddo
762  enddo
763 
764  call read_data (fname, 'PHIS', wk2, no_domain=.true.)
765 ! if(gid==0) call pmaxmin( 'ZS_ncep', wk2, im, jm, 1./grav)
766  do j=js,je
767  do i=is,ie
768  i1 = id1(i,j)
769  i2 = id2(i,j)
770  j1 = jdc(i,j)
771  gz0(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
772  s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
773  enddo
774  enddo
775  call prt_maxmin('ZS_ncep', gz0, is, ie, js, je, 0, 1, 1./grav, master)
776 
777  if ( read_ts ) then ! read skin temperature; could be used for SST
778 
779  call read_data (fname, 'TS', wk2, no_domain=.true.)
780 
781  if ( .not. land_ts ) then
782  allocate ( oro(im,jm) )
783 ! Read NCEP ORO (1; land; 0: ocean; 2: sea_ice)
784  call read_data (fname, 'ORO', oro, no_domain=.true.)
785 
786  do j=1,jm
787  tmean = 0.
788  npt = 0
789  do i=1,im
790  if( abs(oro(i,j)-1.) > 0.5 ) then
791  tmean = tmean + wk2(i,j)
792  npt = npt + 1
793  endif
794  enddo
795 !-------------------------------------------------------
796 ! Replace TS over interior land with zonal mean SST/Ice
797 !-------------------------------------------------------
798  if ( npt /= 0 ) then
799  tmean= tmean / real(npt)
800  do i=1,im
801  if( abs(oro(i,j)-1.) <= 0.5 ) then
802  if ( i==1 ) then
803  i1 = im; i2 = 2
804  elseif ( i==im ) then
805  i1 = im-1; i2 = 1
806  else
807  i1 = i-1; i2 = i+1
808  endif
809  if ( abs(oro(i2,j)-1.)>0.5 ) then ! east side has priority
810  wk2(i,j) = wk2(i2,j)
811  elseif ( abs(oro(i1,j)-1.)>0.5 ) then ! west side
812  wk2(i,j) = wk2(i1,j)
813  else
814  wk2(i,j) = tmean
815  endif
816  endif
817  enddo
818  endif
819  enddo
820  deallocate ( oro )
821  endif ! land_ts
822 
823  if(gid==0) call pmaxmin('SST_ncep', wk2, im, jm, 1.)
824  do j=js,je
825  do i=is,ie
826  i1 = id1(i,j)
827  i2 = id2(i,j)
828  j1 = jdc(i,j)
829  ts(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
830  s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
831  enddo
832  enddo
833  call prt_maxmin('SST_model', ts, is, ie, js, je, 0, 1, 1., master)
834 
835 ! Perform interp to FMS SST format/grid
836  call ncep2fms( wk2 )
837  if(gid==0) call pmaxmin( 'SST_ncep_fms', sst_ncep, i_sst, j_sst, 1.)
838 
839  endif ! read_ts
840 
841  deallocate ( wk2 )
842 
843 ! Read in temperature:
844  allocate ( wk3(im,jm,km) )
845 
846 ! Winds:
847  if ( nudge_winds ) then
848 
849  call read_data (fname, 'U', wk3, no_domain=.true.)
850  if( master ) call pmaxmin( 'U_ncep', wk3, im*jm, km, 1.)
851 
852  do k=1,km
853  do j=js,je
854  do i=is,ie
855  i1 = id1(i,j)
856  i2 = id2(i,j)
857  j1 = jdc(i,j)
858  u(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
859  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
860  enddo
861  enddo
862  enddo
863 
864  call read_data (fname, 'V', wk3, no_domain=.true.)
865  if( master ) call pmaxmin( 'V_ncep', wk3, im*jm, km, 1.)
866  do k=1,km
867  do j=js,je
868  do i=is,ie
869  i1 = id1(i,j)
870  i2 = id2(i,j)
871  j1 = jdc(i,j)
872  v(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
873  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
874  enddo
875  enddo
876  enddo
877 
878  endif
879 
880  if ( nudge_t .or. nudge_virt .or. nudge_q .or. nudge_tpw .or. nudge_hght ) then
881 
882 ! Read in tracers: only sphum at this point
883  call read_data (fname, 'Q', wk3, no_domain=.true.)
884  if(gid==1) call pmaxmin( 'Q_ncep', wk3, im*jm, km, 1.)
885  do k=1,km
886  do j=js,je
887  do i=is,ie
888  i1 = id1(i,j)
889  i2 = id2(i,j)
890  j1 = jdc(i,j)
891  q(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
892  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
893  enddo
894  enddo
895  enddo
896 
897  call read_data (fname, 'T', wk3, no_domain=.true.)
898  if(gid==0) call pmaxmin( 'T_ncep', wk3, im*jm, km, 1.)
899 
900  do k=1,km
901  do j=js,je
902  do i=is,ie
903  i1 = id1(i,j)
904  i2 = id2(i,j)
905  j1 = jdc(i,j)
906  t(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
907  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
908 ! Convert t to virtual temperature:
909  t(i,j,k) = t(i,j,k)*(1.+zvir*q(i,j,k))
910  enddo
911  enddo
912  enddo
913 
914  endif
915 
916  deallocate ( wk3 )
917 
918  nfile = nfile + 1
919 
920  end subroutine get_ncep_analysis
921 
922 
923 
924  subroutine remap_coef
925 
926 ! local:
927  real(FVPRC) :: rdlon(im)
928  real(FVPRC) :: rdlat(jm)
929  real(FVPRC):: a1, b1
930  integer i,j, i1, i2, jc, i0, j0
931 
932  do i=1,im-1
933  rdlon(i) = 1. / (lon(i+1) - lon(i))
934  enddo
935  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
936 
937  do j=1,jm-1
938  rdlat(j) = 1. / (lat(j+1) - lat(j))
939  enddo
940 
941 ! * Interpolate to cubed sphere cell center
942  do 5000 j=js,je
943 
944  do i=is,ie
945 
946  if ( agrid(i,j,1)>lon(im) ) then
947  i1 = im; i2 = 1
948  a1 = (agrid(i,j,1)-lon(im)) * rdlon(im)
949  elseif ( agrid(i,j,1)<lon(1) ) then
950  i1 = im; i2 = 1
951  a1 = (agrid(i,j,1)+2.*pi-lon(im)) * rdlon(im)
952  else
953  do i0=1,im-1
954  if ( agrid(i,j,1)>=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then
955  i1 = i0; i2 = i0+1
956  a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0)
957  go to 111
958  endif
959  enddo
960  endif
961 111 continue
962 
963  if ( agrid(i,j,2)<lat(1) ) then
964  jc = 1
965  b1 = 0.
966  elseif ( agrid(i,j,2)>lat(jm) ) then
967  jc = jm-1
968  b1 = 1.
969  else
970  do j0=1,jm-1
971  if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then
972  jc = j0
973  b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc)
974  go to 222
975  endif
976  enddo
977  endif
978 222 continue
979 
980  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
981  write(*,*) 'gid=', gid, i,j,a1, b1
982  endif
983 
984  s2c(i,j,1) = (1.-a1) * (1.-b1)
985  s2c(i,j,2) = a1 * (1.-b1)
986  s2c(i,j,3) = a1 * b1
987  s2c(i,j,4) = (1.-a1) * b1
988  id1(i,j) = i1
989  id2(i,j) = i2
990  jdc(i,j) = jc
991  enddo !i-loop
992 5000 continue ! j-loop
993 
994  end subroutine remap_coef
995 
996 
997  subroutine ncep2fms( sst )
998  real(FVPRC), intent(in):: sst(im,jm)
999 ! local:
1000  real(FVPRC) :: rdlon(im)
1001  real(FVPRC) :: rdlat(jm)
1002  real(FVPRC):: a1, b1
1003  real(FVPRC):: delx, dely
1004  real(FVPRC):: xc, yc ! "data" location
1005  real(FVPRC):: c1, c2, c3, c4
1006  integer i,j, i1, i2, jc, i0, j0, it, jt
1007 
1008  do i=1,im-1
1009  rdlon(i) = 1. / (lon(i+1) - lon(i))
1010  enddo
1011  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
1012 
1013  do j=1,jm-1
1014  rdlat(j) = 1. / (lat(j+1) - lat(j))
1015  enddo
1016 
1017 ! * Interpolate to "FMS" 1x1 SST data grid
1018 ! lon: 0.5, 1.5, ..., 359.5
1019 ! lat: -89.5, -88.5, ... , 88.5, 89.5
1020 
1021  delx = 360./real(i_sst)
1022  dely = 180./real(j_sst)
1023 
1024  jt = 1
1025  do 5000 j=1,j_sst
1026 
1027  yc = (-90. + dely * (0.5+real(j-1))) * deg2rad
1028  if ( yc<lat(1) ) then
1029  jc = 1
1030  b1 = 0.
1031  elseif ( yc>lat(jm) ) then
1032  jc = jm-1
1033  b1 = 1.
1034  else
1035  do j0=jt,jm-1
1036  if ( yc>=lat(j0) .and. yc<=lat(j0+1) ) then
1037  jc = j0
1038  jt = j0
1039  b1 = (yc-lat(jc)) * rdlat(jc)
1040  go to 222
1041  endif
1042  enddo
1043  endif
1044 222 continue
1045  it = 1
1046 
1047  do i=1,i_sst
1048  xc = delx * (0.5+real(i-1)) * deg2rad
1049  if ( xc>lon(im) ) then
1050  i1 = im; i2 = 1
1051  a1 = (xc-lon(im)) * rdlon(im)
1052  elseif ( xc<lon(1) ) then
1053  i1 = im; i2 = 1
1054  a1 = (xc+2.*pi-lon(im)) * rdlon(im)
1055  else
1056  do i0=it,im-1
1057  if ( xc>=lon(i0) .and. xc<=lon(i0+1) ) then
1058  i1 = i0; i2 = i0+1
1059  it = i0
1060  a1 = (xc-lon(i1)) * rdlon(i0)
1061  go to 111
1062  endif
1063  enddo
1064  endif
1065 111 continue
1066 
1067 ! if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
1068 ! write(*,*) 'gid=', gid, i,j,a1, b1
1069 ! endif
1070  c1 = (1.-a1) * (1.-b1)
1071  c2 = a1 * (1.-b1)
1072  c3 = a1 * b1
1073  c4 = (1.-a1) * b1
1074 ! Interpolated surface pressure
1075  sst_ncep(i,j) = c1*sst(i1,jc ) + c2*sst(i2,jc ) + &
1076  c3*sst(i2,jc+1) + c4*sst(i1,jc+1)
1077  enddo !i-loop
1078 5000 continue ! j-loop
1079 
1080  end subroutine ncep2fms
1081 
1082 
1083  subroutine get_int_hght(h_int, npz, ak, bk, ps, delp, ps0, tv)
1084  integer, intent(in):: npz
1085  real(FVPRC), intent(in):: ak(npz+1), bk(npz+1)
1086  real(FVPRC), intent(in), dimension(is:ie,js:je):: ps, ps0
1087  real(FVPRC), intent(in), dimension(isd:ied,jsd:jed,npz):: delp
1088  real*4, intent(in), dimension(is:ie,js:je,km):: tv
1089  real(FVPRC), intent(out), dimension(is:ie,js:je):: h_int ! g*height
1090 ! local:
1091  real(FVPRC), dimension(is:ie,km+1):: pn0, gz
1092  real(FVPRC):: logp(is:ie)
1093  integer i,j,k
1094 
1095  h_int(:,:) = 1.e25
1096 
1097  do 5000 j=js,je
1098 
1099  do k=1,km+1
1100  do i=is,ie
1101  pn0(i,k) = log( ak0(k) + bk0(k)*ps0(i,j) )
1102  enddo
1103  enddo
1104 !------
1105 ! Model
1106 !------
1107  do i=is,ie
1108  logp(i) = ak(1)
1109  enddo
1110  do k=1,k_trop
1111  do i=is,ie
1112  logp(i) = logp(i) + delp(i,j,k)
1113  enddo
1114  enddo
1115  do i=is,ie
1116  logp(i) = log( logp(i) )
1117  gz(i,km+1) = gz0(i,j) ! Data Surface geopotential
1118  enddo
1119 
1120 ! Linear in log-p interpolation
1121  do i=is,ie
1122  do k=km,1,-1
1123  gz(i,k) = gz(i,k+1) + rdgas*tv(i,j,k)*(pn0(i,k+1)-pn0(i,k))
1124  if ( logp(i)>=pn0(i,k) .and. logp(i)<=pn0(i,k+1) ) then
1125  h_int(i,j) = gz(i,k+1) + (gz(i,k)-gz(i,k+1))*(pn0(i,k+1)-logp(i))/(pn0(i,k+1)-pn0(i,k))
1126  goto 400
1127  endif
1128  enddo
1129 400 continue
1130  enddo
1131 
1132 5000 continue
1133 
1134 
1135  end subroutine get_int_hght
1136 
1137 
1138 
1139  subroutine remap_tq( npz, ak, bk, ps, delp, t, q, &
1140  kmd, ps0, ta, qa, zvir)
1141  integer, intent(in):: npz, kmd
1142  real(FVPRC), intent(in):: zvir
1143  real(FVPRC), intent(in):: ak(npz+1), bk(npz+1)
1144  real(FVPRC), intent(in), dimension(is:ie,js:je):: ps0
1145  real(FVPRC), intent(inout), dimension(is:ie,js:je):: ps
1146  real(FVPRC), intent(in), dimension(isd:ied,jsd:jed,npz):: delp
1147  real*4, intent(in), dimension(is:ie,js:je,kmd):: ta
1148  real*4, intent(in), dimension(is:ie,js:je,kmd):: qa
1149  real*4, intent(out), dimension(is:ie,js:je,npz):: t
1150  real*4, intent(out), dimension(is:ie,js:je,npz):: q
1151 ! local:
1152  real(FVPRC), dimension(is:ie,kmd):: tp, qp
1153  real(FVPRC), dimension(is:ie,kmd+1):: pe0, pn0
1154  real(FVPRC), dimension(is:ie,npz):: qn1
1155  real(FVPRC), dimension(is:ie,npz+1):: pe1, pn1
1156  integer i,j,k
1157 
1158 
1159  do 5000 j=js,je
1160 
1161  do k=1,kmd+1
1162  do i=is,ie
1163  pe0(i,k) = ak0(k) + bk0(k)*ps0(i,j)
1164  pn0(i,k) = log(pe0(i,k))
1165  enddo
1166  enddo
1167 !------
1168 ! Model
1169 !------
1170  do i=is,ie
1171  pe1(i,1) = ak(1)
1172  enddo
1173  do k=1,npz
1174  do i=is,ie
1175  pe1(i,k+1) = pe1(i,k) + delp(i,j,k)
1176  enddo
1177  enddo
1178  do i=is,ie
1179  ps(i,j) = pe1(i,npz+1)
1180  enddo
1181  do k=1,npz+1
1182  do i=is,ie
1183  pn1(i,k) = log(pe1(i,k))
1184  enddo
1185  enddo
1186 
1187  if ( nudge_t .or. nudge_q ) then
1188  do k=1,kmd
1189  do i=is,ie
1190  qp(i,k) = qa(i,j,k)
1191  enddo
1192  enddo
1193  call mappm(kmd, pe0, qp, npz, pe1, qn1, is,ie, 0, kord_data)
1194  do k=1,npz
1195  do i=is,ie
1196  q(i,j,k) = qn1(i,k)
1197  enddo
1198  enddo
1199  endif
1200 
1201  do k=1,kmd
1202  do i=is,ie
1203  tp(i,k) = ta(i,j,k)
1204  enddo
1205  enddo
1206  call mappm(kmd, pn0, tp, npz, pn1, qn1, is,ie, 1, kord_data)
1207 
1208  if ( nudge_t ) then
1209  do k=1,npz
1210  do i=is,ie
1211  t(i,j,k) = qn1(i,k)/(1.+zvir*q(i,j,k))
1212  enddo
1213  enddo
1214  else
1215  do k=1,npz
1216  do i=is,ie
1217  t(i,j,k) = qn1(i,k)
1218  enddo
1219  enddo
1220  endif
1221 
1222 5000 continue
1223 
1224  end subroutine remap_tq
1225 
1226 
1227  subroutine remap_uv(npz, ak, bk, ps, delp, u, v, kmd, ps0, u0, v0)
1228  integer, intent(in):: npz
1229  real(FVPRC), intent(in):: ak(npz+1), bk(npz+1)
1230  real(FVPRC), intent(inout):: ps(is:ie,js:je)
1231  real(FVPRC), intent(in), dimension(isd:ied,jsd:jed,npz):: delp
1232  real*4, intent(inout), dimension(is:ie,js:je,npz):: u, v
1233 !
1234  integer, intent(in):: kmd
1235  real(FVPRC), intent(in):: ps0(is:ie,js:je)
1236  real*4, intent(in), dimension(is:ie,js:je,kmd):: u0, v0
1237 !
1238 ! local:
1239  real(FVPRC), dimension(is:ie,kmd+1):: pe0
1240  real(FVPRC), dimension(is:ie,npz+1):: pe1
1241  real(FVPRC), dimension(is:ie,kmd):: qt
1242  real(FVPRC), dimension(is:ie,npz):: qn1
1243  integer i,j,k
1244 
1245  do 5000 j=js,je
1246 !------
1247 ! Data
1248 !------
1249  do k=1,kmd+1
1250  do i=is,ie
1251  pe0(i,k) = ak0(k) + bk0(k)*ps0(i,j)
1252  enddo
1253  enddo
1254 !------
1255 ! Model
1256 !------
1257  do i=is,ie
1258  pe1(i,1) = ak(1)
1259  enddo
1260  do k=1,npz
1261  do i=is,ie
1262  pe1(i,k+1) = pe1(i,k) + delp(i,j,k)
1263  enddo
1264  enddo
1265  do i=is,ie
1266  ps(i,j) = pe1(i,npz+1)
1267  enddo
1268 !------
1269 ! map u
1270 !------
1271  do k=1,kmd
1272  do i=is,ie
1273  qt(i,k) = u0(i,j,k)
1274  enddo
1275  enddo
1276  call mappm(kmd, pe0, qt, npz, pe1, qn1, is,ie, -1, kord_data)
1277  do k=1,npz
1278  do i=is,ie
1279  u(i,j,k) = qn1(i,k)
1280  enddo
1281  enddo
1282 !------
1283 ! map v
1284 !------
1285  do k=1,kmd
1286  do i=is,ie
1287  qt(i,k) = v0(i,j,k)
1288  enddo
1289  enddo
1290  call mappm(kmd, pe0, qt, npz, pe1, qn1, is,ie, -1, kord_data)
1291  do k=1,npz
1292  do i=is,ie
1293  v(i,j,k) = qn1(i,k)
1294  enddo
1295  enddo
1296 5000 continue
1297 
1298  end subroutine remap_uv
1299 
1300 
1301 
1302  subroutine nwp_nudge_end
1303 
1304  deallocate ( ps_dat )
1305  deallocate ( t_dat )
1306  deallocate ( q_dat )
1307 
1308  if ( nudge_winds ) then
1309  deallocate ( u_dat )
1310  deallocate ( v_dat )
1311  endif
1312 
1313  deallocate ( s2c )
1314  deallocate ( id1 )
1315  deallocate ( id2 )
1316  deallocate ( jdc )
1317 
1318  deallocate ( ak0 )
1319  deallocate ( bk0 )
1320  deallocate ( lat )
1321  deallocate ( lon )
1322 
1323  deallocate ( gz0 )
1324 
1325  end subroutine nwp_nudge_end
1326 
1327 
1328  subroutine get_tc_mask(time, mask)
1329  real(FVPRC) :: slp_mask = 100900. ! crtical SLP to apply mask
1330 ! Input
1331  type(time_type), intent(in):: time
1332  real(FVPRC), intent(inout):: mask(is:ie,js:je)
1333 ! local
1334  real(FVPRC):: pos(2)
1335  real(FVPRC):: slp_o ! sea-level pressure (Pa)
1336  real(FVPRC):: r_vor, p_vor
1337  real(FVPRC):: dist
1338  integer n, i, j
1339 
1340  do 5000 n=1,nstorms ! looop through all storms
1341 !----------------------------------------
1342 ! Obtain slp observation
1343 !----------------------------------------
1344  call get_slp_obs(time, nobs_tc(n), x_obs(1,n), y_obs(1,n), mslp_obs(1,n), mslp_out(i,n), rad_out(1,n), &
1345  time_tc(1,n), pos(1), pos(2), slp_o, r_vor, p_vor)
1346 
1347  if ( slp_o<880.e2 .or. slp_o>min(slp_env,slp_mask) .or. abs(pos(2))*rad2deg>38. ) goto 5000 ! next storm
1348 
1349  if ( r_vor < 30.e3 ) then
1350  r_vor = r_min + (slp_env-slp_o)/20.e2*r_inc ! radius of influence
1351  endif
1352 
1353  do j=js, je
1354  do i=is, ie
1355  dist = great_circle_dist(pos, agrid(i,j,1:2), radius)
1356  if( dist < 5.*r_vor ) then
1357  if ( strong_mask ) then
1358  mask(i,j) = mask(i,j) * ( 1. - exp(-(0.5*dist/r_vor)**2)*min(1.,(slp_env-slp_o)/5.e2) )
1359  else
1360 ! Better analysis data (NCEP 2007 and later) may use a weak mask
1361  mask(i,j) = mask(i,j) * ( 1. - exp(-(0.75*dist/r_vor)**2)*min(1.,(slp_env-slp_o)/10.e2) )
1362  endif
1363  endif
1364  enddo ! i-loop
1365  enddo ! end j-loop
1366 
1367 5000 continue
1368 
1369  end subroutine get_tc_mask
1370 
1371 
1372  subroutine breed_slp_inline(nstep, dt, npz, ak, bk, phis, pe, pk, peln, delp, u, v, pt, q, nwat, zvir)
1373 !------------------------------------------------------------------------------------------
1374 ! Purpose: Vortex-breeding by nudging sea-level-pressure towards single point observations
1375 ! Note: conserve water mass, geopotential, and momentum at the expense of dry air mass
1376 !------------------------------------------------------------------------------------------
1377 ! Input
1378  integer, intent(in):: nstep, npz, nwat
1379  real(FVPRC), intent(in):: dt ! (small) time step in seconds
1380  real(FVPRC), intent(in):: zvir
1381  real(FVPRC), intent(in), dimension(npz+1):: ak, bk
1382  real(FVPRC), intent(in):: phis(isd:ied,jsd:jed)
1383 ! Input/Output
1384  real(FVPRC), intent(inout):: u(isd:ied,jsd:jed+1,npz)
1385  real(FVPRC), intent(inout):: v(isd:ied+1,jsd:jed,npz)
1386  real(FVPRC), intent(inout), dimension(isd:ied,jsd:jed,npz):: delp, pt
1387  real(FVPRC), intent(inout)::q(isd:ied,jsd:jed,npz,*)
1388 
1389  real(FVPRC), intent(inout):: pk(is:ie,js:je, npz+1) ! pe**kappa
1390  real(FVPRC), intent(inout):: pe(is-1:ie+1, npz+1,js-1:je+1) ! edge pressure (pascal)
1391  real(FVPRC), intent(out):: peln(is:ie,npz+1,js:je) ! ln(pe)
1392 ! local
1393  type(time_type):: time
1394  real(FVPRC):: ps(is:ie,js:je)
1395  real(FVPRC):: dist(is:ie,js:je)
1396  real(FVPRC):: tm(is:ie,js:je)
1397  real(FVPRC):: slp(is:ie,js:je)
1398  real(FVPRC):: pos(2)
1399  real(FVPRC):: slp_o ! sea-level pressure (Pa)
1400  real(FVPRC):: p_env
1401  real(FVPRC):: r_vor
1402  real(FVPRC):: relx0, relx, f1, pbreed, pbtop, pkz
1403  real(FVPRC):: p_obs, ratio, p_count, p_sum, mass_sink, delps
1404  real(FVPRC):: p_lo, p_hi, tau_vt
1405  real(FVPRC):: split_time, fac, pdep
1406  integer year, month, day, hour, minute, second
1407  integer n, i, j, k, iq, k0
1408 
1409  if ( nstorms==0 ) then
1410  if(master) write(*,*) 'NO TC data to process'
1411  return
1412  endif
1413 
1414  if ( k_breed==0 ) k_breed = ks
1415 ! k_breed = ks
1416 
1417  k0 = k_breed
1418 
1419 ! Advance (local) time
1420  call get_date(fv_time, year, month, day, hour, minute, second)
1421  if ( year /= year_track_data ) then
1422  if (master) write(*,*) 'Warning: The year in storm track data is not the same as model year'
1423  return
1424  endif
1425  time = fv_time
1426  split_time = calday(year, month, day, hour, minute, second) + dt*real(nstep)/86400.
1427 
1428  do j=js,je
1429 ! ---- Compute ps
1430  do i=is,ie
1431  ps(i,j) = ak(1)
1432  enddo
1433  do k=1,npz
1434  do i=is,ie
1435  ps(i,j) = ps(i,j) + delp(i,j,k)
1436  enddo
1437  enddo
1438 ! Compute lowest layer air temperature:
1439  do i=is,ie
1440  pkz = (pk(i,j,npz+1)-pk(i,j,npz))/(kappa*log(ps(i,j)/(ps(i,j)-delp(i,j,npz))))
1441  tm(i,j) = pkz*pt(i,j,npz)/(cp_air*(1.+zvir*q(i,j,npz,1)))
1442  enddo
1443  enddo
1444 
1445  do k=k_breed+1,npz
1446 
1447  do j=js,je+1
1448  do i=is,ie
1449  u(i,j,k) = u(i,j,k) * (delp(i,j-1,k)+delp(i,j,k))
1450  enddo
1451  enddo
1452  do j=js,je
1453  do i=is,ie+1
1454  v(i,j,k) = v(i,j,k) * (delp(i-1,j,k)+delp(i,j,k))
1455  enddo
1456  enddo
1457  do j=js,je
1458  do i=is,ie
1459  pt(i,j,k) = pt(i,j,k)*(pk(i,j,k+1)-pk(i,j,k))
1460  enddo
1461  enddo
1462 ! Convert tracer moist mixing ratio to mass
1463  do iq=1,nwat
1464  do j=js,je
1465  do i=is,ie
1466  q(i,j,k,iq) = q(i,j,k,iq) * delp(i,j,k)
1467  enddo
1468  enddo
1469  enddo
1470 
1471  enddo
1472 
1473  do 5000 n=1,nstorms ! looop through all storms
1474 
1475 !----------------------------------------
1476 ! Obtain slp observation
1477 !----------------------------------------
1478  call get_slp_obs(time, nobs_tc(n), x_obs(1,n), y_obs(1,n), mslp_obs(1,n), mslp_out(i,n), rad_out(1,n), &
1479  time_tc(1,n), pos(1), pos(2), slp_o, r_vor, p_env, stime=split_time, fact=fac)
1480 
1481  if ( slp_o<87500. .or. slp_o>slp_env .or. abs(pos(2))*rad2deg>40. ) then
1482  goto 5000 ! next storm
1483  endif
1484 
1485  if(nudge_debug .and. master) &
1486  write(*,*) 'Vortex breeding for TC:', n, ' slp=',slp_o/100.,pos(1)*rad2deg,pos(2)*rad2deg
1487 
1488 #ifndef CONST_BREED_DEPTH
1489 ! Determine pbtop (top pressure of vortex breeding)
1490 
1491  if ( slp_o > 1000.e2 ) then
1492  pbtop = 850.e2
1493  else
1494  pbtop = max(125.e2, 850.e2-30.*(1000.e2-slp_o))
1495  endif
1496 
1497  do k=1,npz
1498  pbreed = ak(k) + bk(k)*1000.e2
1499  if ( pbreed>pbtop ) then
1500  k0 = k
1501  exit
1502  endif
1503  enddo
1504  k0 = max(k0, k_breed)
1505 #endif
1506 
1507  do j=js, je
1508  do i=is, ie
1509  dist(i,j) = great_circle_dist( pos, agrid(i,j,1:2), radius)
1510  enddo
1511  enddo
1512 ! ---- Compute slp
1513  do j=js,je
1514  do i=is,ie
1515  slp(i,j) = ps(i,j)*exp(phis(i,j)/(rdgas*(tm(i,j)+3.25e-3*phis(i,j)/grav)))
1516  enddo
1517  enddo
1518 
1519 
1520  if ( r_vor < 30.e3 .or. p_env<900.e2 ) then
1521 
1522 ! Compute r_vor & p_env
1523  r_vor = r_min + (slp_env-slp_o)/25.e2*r_inc
1524 
1525 123 continue
1526  p_count = 0.
1527  p_sum = 0.
1528  do j=js, je
1529  do i=is, ie
1530  if( dist(i,j)<(r_vor+del_r) .and. dist(i,j)>r_vor .and. phis(i,j)<200.*grav ) then
1531  p_count = p_count + 1.
1532  p_sum = p_sum + slp(i,j)
1533  endif
1534  enddo
1535  enddo
1536 
1537  call mp_reduce_sum(p_count)
1538 
1539  if ( p_count<32. ) then
1540  if(nudge_debug .and. master) write(*,*) p_count, 'Skipping obs: too few p_count'
1541  goto 5000
1542  endif
1543 
1544  call mp_reduce_sum(p_sum)
1545  p_env = p_sum / p_count
1546 
1547  if(nudge_debug .and. master) write(*,*) 'Environmental SLP=', p_env/100., ' computed radius=', r_vor/1.e3
1548 
1549  if ( p_env>1020.e2 .or. p_env<900.e2 ) then
1550  if( nudge_debug ) then
1551  if(master) write(*,*) 'Environmental SLP out of bound; skipping obs. p_count=', p_count, p_sum
1552  call prt_maxmin('SLP_breeding', slp, is, ie, js, je, 0, 1, 0.01, master)
1553  endif
1554  goto 5000
1555  endif
1556 
1557  endif
1558 
1559 ! if ( p_env < max(slp_o, 1000.E2) ) then
1560  if ( p_env < max(1000.e2, slp_o+250.0) ) then
1561  if(nudge_debug .and. master) then
1562  write(*,*) 'Computed environmental SLP too low'
1563  write(*,*) ' ', p_env/100., slp_o/100.,pos(1)*rad2deg, pos(2)*rad2deg
1564  endif
1565  if ( r_vor < 750.e3 ) then
1566  r_vor = r_vor + del_r
1567  if(nudge_debug .and. master) write(*,*) 'Vortex radius (km) increased to:', r_vor/1.e3
1568  goto 123
1569  else
1570  p_env = max( slp_o + 250.0, 1000.e2)
1571  endif
1572  endif
1573 
1574 ! make tau linear function of SLP_OBS
1575 ! Stronger nudging for weak cyclones
1576  tau_vt = tau_vortex + 6.*(980.e2-slp_o)/100.
1577  tau_vt = max(dt, tau_vt)
1578 
1579  if ( time_track ) then
1580  relx0 = min(1., fac*dt/tau_vt)
1581  else
1582  relx0 = min(1., dt/tau_vt)
1583  endif
1584  mass_sink = 0.
1585  do j=js, je
1586  do i=is, ie
1587  if( dist(i,j) < r_vor .and. phis(i,j)<200.*grav ) then
1588  f1 = dist(i,j)/r_vor
1589  relx = relx0*exp( -4.*f1**2 )
1590 ! Compute p_obs: assuming local radial distributions of slp are Gaussian
1591  p_hi = p_env - (p_env-slp_o) * exp( -5.0*f1**2 ) ! upper bound
1592  p_lo = p_env - (p_env-slp_o) * exp( -2.0*f1**2 ) ! lower bound
1593 
1594  if ( ps(i,j) > p_hi ) then
1595 ! Under-development:
1596  delps = relx*(ps(i,j) - p_hi) ! Note: ps is used here to prevent
1597  ! over deepening over terrain
1598  elseif ( slp(i,j) < p_lo ) then
1599 ! Over-development:
1600  delps = relx*(slp(i,j) - p_lo) ! Note: slp is used here
1601  else
1602  goto 400 ! do nothing; proceed to next storm
1603  endif
1604 
1605  mass_sink = mass_sink + delps*area(i,j)
1606 
1607 !==========================================================================================
1608  if ( delps > 0. ) then
1609  pbreed = ak(1)
1610  do k=1,k0
1611  pbreed = pbreed + delp(i,j,k)
1612  enddo
1613  f1 = 1. - delps/(ps(i,j)-pbreed)
1614  do k=k0+1,npz
1615  delp(i,j,k) = delp(i,j,k)*f1
1616  enddo
1617  else
1618  do k=npz,k0+2,-1
1619  if ( abs(delps) < 1. ) then
1620  delp(i,j,k) = delp(i,j,k) - delps
1621  go to 400
1622  else
1623 ! pdep = max(1.0, min(abs(0.5*delps), 0.020*delp(i,j,k)))
1624  pdep = max(1.0, min(abs(0.4*delps), 0.025*delp(i,j,k)))
1625  pdep = sign( min(pdep, abs(delps)), delps )
1626  delp(i,j,k) = delp(i,j,k) - pdep
1627  delps = delps - pdep
1628  endif
1629  enddo
1630 ! Final stage: put all remaining correction term to the (npz-k_breed) layer
1631  delp(i,j,k0+1) = delp(i,j,k0+1) - delps
1632  endif
1633 !==========================================================================================
1634 
1635  endif
1636 400 continue
1637  enddo ! end i-loop
1638  enddo ! end j-loop
1639 
1640  call mp_reduce_sum(mass_sink)
1641  if ( abs(mass_sink)<1.e-40 ) goto 5000
1642 
1643  p_sum = 0.
1644  do j=js, je
1645  do i=is, ie
1646  if( dist(i,j)<6.*r_vor .and. dist(i,j)>r_vor+del_r ) then
1647  p_sum = p_sum + area(i,j)
1648  endif
1649  enddo
1650  enddo
1651 
1652  call mp_reduce_sum(p_sum)
1653  mass_sink = mass_sink / p_sum ! mean delta pressure to be added back to the environment to conserve mass
1654  if(master .and. nudge_debug) write(*,*) 'TC#',n, 'Mass tele-ported (pa)=', mass_sink
1655 
1656  do j=js, je
1657  do i=is, ie
1658  if( dist(i,j)<6.*r_vor .and. dist(i,j)>r_vor+del_r ) then
1659  pbreed = ak(1)
1660  do k=1,k_breed
1661  pbreed = pbreed + delp(i,j,k)
1662  enddo
1663  f1 = 1. + mass_sink/(ps(i,j)-pbreed)
1664  do k=k_breed+1,npz
1665  delp(i,j,k) = delp(i,j,k)*f1
1666  enddo
1667  endif
1668  enddo
1669  enddo
1670 
1671 ! ---- re-compute ps
1672  do j=js,je
1673  do i=is,ie
1674  ps(i,j) = ak(1)
1675  enddo
1676  do k=1,npz
1677  do i=is,ie
1678  ps(i,j) = ps(i,j) + delp(i,j,k)
1679  enddo
1680  enddo
1681  enddo
1682 
1683 5000 continue
1684 
1685  call mpp_update_domains(delp, domain, complete=.true.)
1686 
1687  do j=js-1,je+1
1688  do i=is-1,ie+1
1689  pe(i,1,j) = ak(1)
1690  enddo
1691  do k=2,npz+1
1692  do i=is-1,ie+1
1693  pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
1694  enddo
1695  enddo
1696  enddo
1697 
1698  do k=k_breed+1,npz+1
1699  do j=js,je
1700  do i=is,ie
1701  peln(i,k,j) = log(pe(i,k,j))
1702  pk(i,j,k) = pe(i,k,j)**kappa
1703  enddo
1704  enddo
1705  enddo
1706 
1707 
1708  do k=k_breed+1,npz
1709  do j=js,je+1
1710  do i=is,ie
1711  u(i,j,k) = u(i,j,k) / (delp(i,j-1,k)+delp(i,j,k))
1712  enddo
1713  enddo
1714  do j=js,je
1715  do i=is,ie+1
1716  v(i,j,k) = v(i,j,k) / (delp(i-1,j,k)+delp(i,j,k))
1717  enddo
1718  enddo
1719  do j=js,je
1720  do i=is,ie
1721  pt(i,j,k) = pt(i,j,k) / (pk(i,j,k+1)-pk(i,j,k))
1722  enddo
1723  enddo
1724  enddo
1725 
1726 ! Convert tracer mass back to moist mixing ratio
1727  do iq=1,nwat
1728  do k=k_breed+1,npz
1729  do j=js,je
1730  do i=is,ie
1731  q(i,j,k,iq) = q(i,j,k,iq) / delp(i,j,k)
1732  enddo
1733  enddo
1734  enddo
1735  enddo
1736 
1737  call mpp_update_domains(pt, domain, complete=.true.)
1738 
1739  end subroutine breed_slp_inline
1740 
1741 
1742  subroutine breed_slp(time, dt, npz, ak, bk, ps, phis, delp, ua, va, u_dt, v_dt, pt, q, nwat, zvir)
1743 !------------------------------------------------------------------------------------------
1744 ! Purpose: Vortex-breeding by nudging sea-level-pressure towards single point observations
1745 ! Note: conserve water mass, geopotential, and momentum at the expense of dry air mass
1746 !------------------------------------------------------------------------------------------
1747 ! Input
1748  type(time_type), intent(in):: time
1749  integer, intent(in):: npz, nwat
1750  real(FVPRC), intent(in):: dt ! time step in seconds
1751  real(FVPRC), intent(in):: zvir
1752  real(FVPRC), intent(in), dimension(npz+1):: ak, bk
1753  real(FVPRC), intent(in):: phis(isd:ied,jsd:jed)
1754  real(FVPRC), intent(in):: ps(isd:ied,jsd:jed)
1755 ! Input/Output
1756  real(FVPRC), intent(inout), dimension(isd:ied,jsd:jed,npz):: delp, pt, ua, va, u_dt, v_dt
1757  real(FVPRC), intent(inout)::q(isd:ied,jsd:jed,npz,nwat)
1758 ! local
1759  real(FVPRC):: dist(is:ie,js:je)
1760  real(FVPRC):: slp(is:ie,js:je)
1761  real(FVPRC):: p0(npz+1), p1(npz+1), dm(npz), tvir(npz)
1762  real(FVPRC):: pos(2)
1763  real(FVPRC):: slp_o ! sea-level pressure (Pa)
1764  real(FVPRC):: p_env
1765  real(FVPRC):: r_vor
1766  real(FVPRC):: relx0, relx, f1, rdt
1767  real(FVPRC):: p_obs, ratio, p_count, p_sum, mass_sink, delps
1768  real(FVPRC):: p_lo, p_hi
1769  integer n, i, j, k, iq
1770 
1771  if ( nstorms==0 ) then
1772  if(master) write(*,*) 'NO TC data to process'
1773  return
1774  endif
1775 
1776  rdt = 1./dt
1777  relx0 = min(1., dt/tau_vortex)
1778 
1779  do j=js, je
1780  do i=is, ie
1781  slp(i,j) = ps(i,j)*exp(phis(i,j)/(rdgas*(pt(i,j,npz)+3.25e-3*phis(i,j)/grav)))
1782  enddo
1783  enddo
1784 
1785  do 5000 n=1,nstorms ! looop through all storms
1786 
1787 !----------------------------------------
1788 ! Obtain slp observation
1789 !----------------------------------------
1790  call get_slp_obs(time, nobs_tc(n), x_obs(1,n), y_obs(1,n), mslp_obs(1,n), mslp_out(i,n), rad_out(1,n), &
1791  time_tc(1,n), pos(1), pos(2), slp_o, r_vor, p_env)
1792 
1793  if ( slp_o<87000. .or. slp_o>slp_env .or. abs(pos(2))*rad2deg>40. ) then
1794  goto 5000 ! next storm
1795  endif
1796 
1797  if(nudge_debug .and. master) &
1798  write(*,*) 'Vortex breeding for TC:', n, ' slp=',slp_o/100.,pos(1)*rad2deg,pos(2)*rad2deg
1799 
1800  do j=js, je
1801  do i=is, ie
1802  dist(i,j) = great_circle_dist( pos, agrid(i,j,1:2), radius)
1803  enddo
1804  enddo
1805 
1806  if ( r_vor < 30.e3 .or. p_env<900.e2 ) then
1807 
1808 ! Compute r_vor & p_env
1809  r_vor = r_min + (slp_env-slp_o)/25.e2*r_inc
1810 
1811 123 continue
1812  p_count = 0.
1813  p_sum = 0.
1814  do j=js, je
1815  do i=is, ie
1816  if( dist(i,j)<(r_vor+del_r) .and. dist(i,j)>r_vor .and. phis(i,j)<250.*grav ) then
1817  p_count = p_count + 1.
1818  p_sum = p_sum + slp(i,j)
1819  endif
1820  enddo
1821  enddo
1822 
1823  call mp_reduce_sum(p_count)
1824 
1825  if ( p_count<32. ) then
1826  if(nudge_debug .and. master) write(*,*) p_count, 'Skipping obs: too few p_count'
1827  goto 5000
1828  endif
1829 
1830  call mp_reduce_sum(p_sum)
1831  p_env = p_sum / p_count
1832 
1833  if(nudge_debug .and. master) write(*,*) 'Environmental SLP=', p_env/100., ' computed radius=', r_vor/1.e3
1834 
1835  if ( p_env>1020.e2 .or. p_env<900.e2 ) then
1836  if( nudge_debug ) then
1837  if(master) write(*,*) 'Environmental SLP out of bound; skipping obs. p_count=', p_count, p_sum
1838  call prt_maxmin('SLP_breeding', slp, is, ie, js, je, 0, 1, 0.01, master)
1839  endif
1840  goto 5000
1841  endif
1842 
1843  endif
1844 
1845  if ( p_env < max(slp_o, 1000.e2) ) then
1846  if(nudge_debug .and. master) then
1847  write(*,*) 'Computed environmental SLP too low'
1848  write(*,*) ' ', p_env/100., slp_o/100.,pos(1)*rad2deg, pos(2)*rad2deg
1849 ! This is rare. But it does happen in W. Pacific (e.g., early development stage, Talim 2005); make radius larger
1850  endif
1851  if ( r_vor < 500.e3 ) then
1852  r_vor = r_vor + del_r
1853  if(nudge_debug .and. master) write(*,*) 'Vortex radius (km) increased to:', r_vor/1.e3
1854  goto 123
1855  else
1856  p_env = max(slp_o + 200., 1000.e2)
1857  endif
1858  endif
1859 
1860  mass_sink = 0.
1861  do j=js, je
1862  do i=is, ie
1863  if( dist(i,j) < r_vor .and. phis(i,j)<250.*grav ) then
1864  p0(ks+1) = ak(ks+1)
1865  do k=ks+1,npz
1866  tvir(k) = pt(i,j,k) * (1.+ zvir*q(i,j,k,1))
1867  p0(k+1) = p0(k) + delp(i,j,k)
1868  enddo
1869 !===============================================================================================
1870  f1 = dist(i,j)/r_vor
1871 ! Compute p_obs: assuming local radial distributions of slp are Gaussian
1872  p_hi = p_env - (p_env-slp_o) * exp( -5.0*f1**2 ) ! upper bound
1873  p_lo = p_env - (p_env-slp_o) * exp( -2.0*f1**2 ) ! lower bound
1874  relx = relx0 * exp( -4.*f1**2 )
1875 ! Compute p_obs: assuming local radial distributions of slp are Gaussian
1876 
1877  if ( ps(i,j) > p_hi ) then
1878 ! under-development
1879  delps = relx*(ps(i,j) - p_hi) ! Note: ps is used here to prevent
1880  ! over deepening over terrain
1881  elseif ( slp(i,j) < p_lo ) then
1882 ! over-development
1883  delps = relx*(slp(i,j) - p_lo) ! Note: slp is used here
1884 
1885 ! if ( slp(i,j) < slp_o ) then
1886 ! delps = min(delps, 0.333*(1.+2.*f1**2)*(slp(i,j)-slp_o))
1887 ! endif
1888  else
1889 ! Leave the model alone If the ps/slp is in between [p_lo,p_hi]
1890  goto 400 ! do nothing; proceed to next storm
1891  endif
1892 !===============================================================================================
1893 
1894  mass_sink = mass_sink + delps*area(i,j)
1895  p1(ks+1) = p0(ks+1)
1896  do k=ks+1,npz
1897  dm(k) = delp(i,j,k) - delps*(bk(k+1)-bk(k))
1898  p1(k+1) = p1(k) + dm(k)
1899  ratio = delp(i,j,k)/dm(k)
1900  delp(i,j,k) = dm(k)
1901  do iq=1,nwat
1902  q(i,j,k,iq) = q(i,j,k,iq) * ratio
1903  enddo
1904 ! Recompute pt by preserving height
1905  pt(i,j,k) = tvir(k)*log(p0(k+1)/p0(k))/(log(p1(k+1)/p1(k))*(1.+zvir*q(i,j,k,1)))
1906 ! Momentum conserving re-scaling
1907  u_dt(i,j,k) = u_dt(i,j,k) + ua(i,j,k)*(ratio - 1.)*rdt
1908  v_dt(i,j,k) = v_dt(i,j,k) + va(i,j,k)*(ratio - 1.)*rdt
1909  ua(i,j,k) = ua(i,j,k) * ratio
1910  va(i,j,k) = va(i,j,k) * ratio
1911  enddo
1912  endif
1913 400 continue
1914  enddo ! end i-loop
1915  enddo ! end j-loop
1916 
1917  call mp_reduce_sum(mass_sink)
1918  if ( mass_sink==0. ) goto 5000
1919 
1920  p_sum = 0.
1921  do j=js, je
1922  do i=is, ie
1923  if( dist(i,j)<(6.*r_vor+del_r) .and. dist(i,j)>r_vor+del_r ) then
1924  p_sum = p_sum + area(i,j)
1925  endif
1926  enddo
1927  enddo
1928 
1929  call mp_reduce_sum(p_sum)
1930  mass_sink = mass_sink / p_sum ! mean delta pressure to be added back to the environment to conserve mass
1931  if(master .and. nudge_debug) write(*,*) 'TC#',n, 'Mass tele-ported (pa)=', mass_sink
1932 
1933  do j=js, je
1934  do i=is, ie
1935  if( dist(i,j)<(6.*r_vor+del_r) .and. dist(i,j)>r_vor+del_r ) then
1936  p0(ks+1) = ak(ks+1)
1937  do k=ks+1,npz
1938  tvir(k) = pt(i,j,k) * (1.+ zvir*q(i,j,k,1))
1939  p0(k+1) = p0(k) + delp(i,j,k)
1940  enddo
1941  p1(ks+1) = p0(ks+1)
1942  do k=ks+1,npz
1943  dm(k) = delp(i,j,k) + (bk(k+1)-bk(k))*mass_sink
1944  p1(k+1) = p1(k) + dm(k)
1945  ratio = delp(i,j,k)/dm(k)
1946  delp(i,j,k) = dm(k)
1947  do iq=1,nwat
1948  q(i,j,k,iq) = q(i,j,k,iq) * ratio
1949  enddo
1950  pt(i,j,k) = tvir(k)*log(p0(k+1)/p0(k))/(log(p1(k+1)/p1(k))*(1.+zvir*q(i,j,k,1)))
1951  u_dt(i,j,k) = u_dt(i,j,k) + ua(i,j,k)*(ratio - 1.)*rdt
1952  v_dt(i,j,k) = v_dt(i,j,k) + va(i,j,k)*(ratio - 1.)*rdt
1953  ua(i,j,k) = ua(i,j,k) * ratio
1954  va(i,j,k) = va(i,j,k) * ratio
1955  enddo
1956  endif
1957  enddo
1958  enddo
1959 
1960 5000 continue
1961 
1962  end subroutine breed_slp
1963 
1964 
1965  subroutine get_slp_obs(time, nobs, lon_obs, lat_obs, mslp, slp_out, r_out, time_obs, &
1966  x_o, y_o, slp_o, r_vor, p_vor, stime, fact)
1967 ! Input
1968  type(time_type), intent(in):: time
1969  integer, intent(in):: nobs ! number of observations in this particular storm
1970  real*4, intent(in):: lon_obs(nobs)
1971  real*4, intent(in):: lat_obs(nobs)
1972  real*4, intent(in):: mslp(nobs) ! observed SLP in pa
1973  real*4, intent(in):: slp_out(nobs) ! slp at r_out
1974  real*4, intent(in):: r_out(nobs) !
1975  real*4, intent(in):: time_obs(nobs)
1976  real(FVPRC), optional, intent(in):: stime
1977  real(FVPRC), optional, intent(out):: fact
1978 ! Output
1979  real(FVPRC), intent(out):: x_o , y_o ! position of the storm center
1980  real(FVPRC), intent(out):: slp_o ! Observed sea-level-pressure (pa)
1981  real(FVPRC), intent(out):: r_vor, p_vor
1982 ! Internal:
1983  real(FVPRC):: p1(2), p2(2)
1984  real(FVPRC) time_model
1985  real(FVPRC) fac
1986  integer year, month, day, hour, minute, second, n
1987 
1988  slp_o = -100000.
1989  x_o = -100.*pi
1990  y_o = -100.*pi
1991  p_vor = -1.e10
1992  r_vor = -1.e10
1993 
1994  if ( present(stime) ) then
1995  time_model = stime
1996  else
1997  call get_date(time, year, month, day, hour, minute, second)
1998 
1999  if ( year /= year_track_data ) then
2000  if (master) write(*,*) 'Warning: The year in storm track data is not the same as model year'
2001  return
2002  endif
2003 
2004  time_model = calday(year, month, day, hour, minute, second)
2005 ! if(nudge_debug .and. master) write(*,*) 'Model:', time_model, year, month, day, hour, minute, second
2006  endif
2007 
2008  if ( time_model <= time_obs(1) .or. time_model >= time_obs(nobs) ) then
2009 ! if(nudge_debug .and. master) write(*,*) ' No valid TC data'
2010  return
2011  else
2012  do n=1,nobs-1
2013  if( time_model >= time_obs(n) .and. time_model <= time_obs(n+1) ) then
2014  fac = (time_model-time_obs(n)) / (time_obs(n+1)-time_obs(n))
2015  slp_o = mslp(n) + ( mslp(n+1)- mslp(n)) * fac
2016 ! Trajectory interpolation:
2017 #ifdef LINEAR_TRAJ
2018 ! Linear in (lon,lat) space
2019  x_o = lon_obs(n) + (lon_obs(n+1)-lon_obs(n)) * fac
2020  y_o = lat_obs(n) + (lat_obs(n+1)-lat_obs(n)) * fac
2021 #else
2022  p1(1) = lon_obs(n); p1(2) = lat_obs(n)
2023  p2(1) = lon_obs(n+1); p2(2) = lat_obs(n+1)
2024  call intp_great_circle(fac, p1, p2, x_o, y_o)
2025 #endif
2026 !----------------------------------------------------------------------
2027  if ( present(fact) ) fact = 1. + 0.5*cos(fac*2.*pi)
2028 ! Additional data from the extended best track
2029 ! if ( slp_out(n)>0. .and. slp_out(n+1)>0. .and. r_out(n)>0. .and. r_out(n+1)>0. ) then
2030 ! p_vor = slp_out(n) + ( slp_out(n+1) - slp_out(n)) * fac
2031 ! r_vor = r_out(n) + ( r_out(n+1) - r_out(n)) * fac
2032 ! endif
2033  return
2034  endif
2035  enddo
2036  endif
2037 
2038  end subroutine get_slp_obs
2039 
2040 
2041  subroutine slp_obs_init
2042  integer:: unit, n, nobs
2043  character(len=3):: GMT
2044  character(len=9):: ts_name
2045  character(len=19):: comment
2046  integer:: mmddhh, yr, year, month, day, hour, MPH, islp
2047  integer:: it, i1, i2, p_ring, d_ring
2048  real(FVPRC):: lon_deg, lat_deg, cald, slp, mps
2049 
2050  nobs_tc(:) = 0
2051  time_tc(:,:) = 0.
2052  mslp_obs(:,:) = -100000.
2053  x_obs(:,:) = - 100.*pi
2054  y_obs(:,:) = - 100.*pi
2055 
2056  mslp_out(:,:) = -1.e10
2057  rad_out(:,:) = -1.e10
2058 
2059  if( track_file_name == "No_File_specified" ) then
2060  if(master) write(*,*) 'No TC track file specified'
2061  return
2062  else
2063  unit = 98
2064  open( unit, file=track_file_name)
2065  endif
2066 
2067  read(unit, *) year
2068  if(master) write(*,*) 'Reading TC track data for YEAR=', year
2069 
2070  year_track_data = year
2071 
2072  nstorms = 0
2073  nobs = 0
2074  month = 99
2075 
2076  if ( ibtrack ) then
2077 
2078 !---------------------------------------------------------------
2079 ! The data format is from Ming Zhoa's processed ibTrack datasets
2080 !---------------------------------------------------------------
2081 
2082  read(unit, *) ts_name, nobs, yr, month, day, hour
2083 
2084  if ( yr /= year ) then
2085  if(master) write(*, *) 'Year inconsistency found !!!'
2086  call mpp_error(fatal,'==> Error in reading best track data')
2087  endif
2088 
2089  do while ( ts_name=='start' )
2090 
2091  nstorms = nstorms + 1
2092  nobs_tc(nstorms) = nobs ! observation count for this storm
2093  if(master) write(*, *) 'Read Data for TC#', nstorms, nobs
2094 
2095  do it=1, nobs
2096  read(unit, *) lon_deg, lat_deg, mps, slp, yr, month, day, hour
2097 ! if ( yr /= year ) then
2098 ! if(master) write(*, *) 'Extended to year + 1', yr
2099 ! endif
2100  cald = calday(yr, month, day, hour, 0, 0)
2101  time_tc(it,nstorms) = cald
2102  if(master) write(*, 100) cald, month, day, hour, lon_deg, lat_deg, mps, slp
2103 
2104  mslp_obs(it,nstorms) = 100.*slp
2105  y_obs(it,nstorms) = lat_deg * deg2rad
2106  x_obs(it,nstorms) = lon_deg * deg2rad
2107  enddo
2108 
2109  read(unit, *) ts_name, nobs, yr, month, day, hour
2110  enddo
2111 100 format(1x, f9.2, 1x, i3, 1x, i2, 1x, i2, 1x, f6.1, 1x, f6.1, 1x, f4.1, 1x, f6.1)
2112 
2113  else
2114 
2115  do while ( month /= 0 )
2116 
2117  read(unit, *) month, day, hour, gmt, lat_deg, lon_deg, mph, islp, comment
2118 
2119  select case (month)
2120 
2121  case (99) ! New storm record to start
2122  nstorms = nstorms + 1
2123  nobs = 0
2124  if(master) write(*, *) 'Reading data for TC#', nstorms, comment
2125  case ( 0) ! end of record
2126  if(master) write(*, *) 'End of record reached'
2127  case default
2128  nobs = nobs + 1
2129  cald = calday(year, month, day, hour, 0, 0)
2130  time_tc(nobs,nstorms) = cald
2131  nobs_tc(nstorms) = nobs ! observation count for this storm
2132 
2133  if(master) write(*, 200) nobs, cald, month, day, hour, gmt, lat_deg, lon_deg, mph, islp, comment
2134  mslp_obs(nobs,nstorms) = 100. * real(islp)
2135  y_obs(nobs,nstorms) = lat_deg * deg2rad
2136  if ( gmt == 'GMT' ) then
2137 ! Transfrom x from (-180 , 180) to (0, 360) then to radian
2138  if ( lon_deg < 0 ) then
2139  x_obs(nobs,nstorms) = (360.+lon_deg) * deg2rad
2140  else
2141  x_obs(nobs,nstorms) = (360.-lon_deg) * deg2rad
2142  endif
2143  elseif ( gmt == 'PAC' ) then ! Pacific storms
2144  x_obs(nobs,nstorms) = lon_deg * deg2rad
2145  endif
2146  end select
2147 
2148  enddo
2149 
2150  endif
2151 
2152  close(unit)
2153 
2154  if(master) then
2155  write(*,*) 'TC vortex breeding: total storms=', nstorms
2156  if ( nstorms/=0 ) then
2157  do n=1,nstorms
2158  write(*, *) 'TC#',n, ' contains ', nobs_tc(n),' observations'
2159  enddo
2160  endif
2161  endif
2162 
2163 200 format(i3, 1x,f9.4, 1x, i2, 1x, i2, 1x, i2, 1x, a3, f5.1, 1x, f5.1, 1x, i3, 1x, i4, 1x, a19)
2164 
2165  end subroutine slp_obs_init
2166 
2167 
2168  real(FVPRC) function calday(year, month, day, hour, minute, sec)
2169 ! For time interpolation; Julian day (0 to 365 for non-leap year)
2170 ! input:
2171  integer, intent(in):: year, month, day, hour
2172  integer, intent(in):: minute, sec
2173 ! Local:
2174  integer n, m, ds, nday
2175  real(FVPRC) tsec
2176  integer days(12)
2177  data days /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
2178 
2179  ds = day - 1
2180 
2181  if( month /= 1 ) then
2182  do m=1, month-1
2183  if( m==2 .and. leap_year(year) ) then
2184  ds = ds + 29
2185  else
2186  ds = ds + days(m)
2187  endif
2188  enddo
2189  endif
2190 
2191  if ( leap_year(year_track_data) ) then
2192  nday = 366
2193  else
2194  nday = 365
2195  endif
2196 
2197  calday = real((year-year_track_data)*nday + ds) + real(hour*3600 + minute*60 + sec)/86400.
2198 
2199  end function calday
2200 
2201 
2202  logical function leap_year(ny)
2203  integer, intent(in):: ny
2204 !
2205 ! Determine if year ny is a leap year
2206 ! Author: S.-J. Lin
2207  integer ny00
2208 !
2209 ! No leap years prior to 0000
2210 !
2211  parameter( ny00 = 0000 ) ! The threshold for starting leap-year
2212 
2213  if( ny >= ny00 ) then
2214  if( mod(ny,100) == 0. .and. mod(ny,400) == 0. ) then
2215  leap_year = .true.
2216  elseif( mod(ny,4) == 0. .and. mod(ny,100) /= 0. ) then
2217  leap_year = .true.
2218  else
2219  leap_year = .false.
2220  endif
2221  else
2222  leap_year = .false.
2223  endif
2224 
2225  end function leap_year
2226 
2227 
2228  subroutine pmaxmin( qname, a, imax, jmax, fac )
2229 
2230  character*(*) qname
2231  integer imax, jmax
2232  integer i, j
2233  real(FVPRC) a(imax,jmax)
2234 
2235  real(FVPRC) qmin(jmax), qmax(jmax)
2236  real(FVPRC) pmax, pmin
2237  real(FVPRC) fac ! multiplication factor
2238 
2239  do j=1,jmax
2240  pmax = a(1,j)
2241  pmin = a(1,j)
2242  do i=2,imax
2243  pmax = max(pmax, a(i,j))
2244  pmin = min(pmin, a(i,j))
2245  enddo
2246  qmax(j) = pmax
2247  qmin(j) = pmin
2248  enddo
2249 !
2250 ! Now find max/min of amax/amin
2251 !
2252  pmax = qmax(1)
2253  pmin = qmin(1)
2254  do j=2,jmax
2255  pmax = max(pmax, qmax(j))
2256  pmin = min(pmin, qmin(j))
2257  enddo
2258 
2259  write(*,*) qname, ' max = ', pmax*fac, ' min = ', pmin*fac
2260 
2261  end subroutine pmaxmin
2262 
2263 
2264  subroutine del2_uv(du, dv, cd, kmd, ntimes)
2265 ! This routine is for filtering the wind tendency
2266  integer, intent(in):: kmd
2267  integer, intent(in):: ntimes
2268  real(FVPRC), intent(in):: cd ! cd = K * da_min; 0 < K < 0.25
2269  real(FVPRC), intent(inout):: du(is:ie,js:je,kmd)
2270  real(FVPRC), intent(inout):: dv(is:ie,js:je,kmd)
2271 ! local:
2272  real(FVPRC), dimension(is:ie,js:je,kmd):: v1, v2, v3
2273  integer i,j,k
2274 
2275 ! transform to 3D Cartesian:
2276  do k=1,kmd
2277  do j=js,je
2278  do i=is,ie
2279  v1(i,j,k) = du(i,j,k)*vlon(i,j,1) + dv(i,j,k)*vlat(i,j,1)
2280  v2(i,j,k) = du(i,j,k)*vlon(i,j,2) + dv(i,j,k)*vlat(i,j,2)
2281  v3(i,j,k) = du(i,j,k)*vlon(i,j,3) + dv(i,j,k)*vlat(i,j,3)
2282  enddo
2283  enddo
2284  enddo
2285 
2286 ! Filter individual components as scalar:
2287  call del2_scalar( v1(is,js,1), cd, kmd, ntimes )
2288  call del2_scalar( v2(is,js,1), cd, kmd, ntimes )
2289  call del2_scalar( v3(is,js,1), cd, kmd, ntimes )
2290 
2291 ! Convert back to lat-lon components:
2292  do k=1,kmd
2293  do j=js,je
2294  do i=is,ie
2295  du(i,j,k) = v1(i,j,k)*vlon(i,j,1) + v2(i,j,k)*vlon(i,j,2) + v3(i,j,k)*vlon(i,j,3)
2296  dv(i,j,k) = v1(i,j,k)*vlat(i,j,1) + v2(i,j,k)*vlat(i,j,2) + v3(i,j,k)*vlat(i,j,3)
2297  enddo
2298  enddo
2299  enddo
2300 
2301  end subroutine del2_uv
2302 
2303  subroutine del2_scalar(qdt, cd, kmd, ntimes)
2304 ! This routine is for filtering the physics tendency
2305  integer, intent(in):: kmd
2306  integer, intent(in):: ntimes
2307  real(FVPRC), intent(in):: cd ! cd = K * da_min; 0 < K < 0.25
2308  real(FVPRC), intent(inout):: qdt(is:ie,js:je,kmd)
2309 ! local:
2310  real(FVPRC):: q(isd:ied,jsd:jed,kmd)
2311  real(FVPRC):: fx(isd:ied+1,jsd:jed), fy(isd:ied,jsd:jed+1)
2312  integer i,j,k, n, nt
2313  real(FVPRC) :: damp
2314 
2315  damp = cd * da_min
2316 
2317  do k=1,kmd
2318  do j=js,je
2319  do i=is,ie
2320  q(i,j,k) = qdt(i,j,k)
2321  enddo
2322  enddo
2323  enddo
2324  call timing_on('COMM_TOTAL')
2325  call mpp_update_domains(q, domain, complete=.true.)
2326  call timing_off('COMM_TOTAL')
2327 
2328  do n=1,ntimes
2329 
2330  nt = ntimes - n
2331 
2332  do k=1,kmd
2333 
2334  if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 1)
2335  do j=js-nt,je+nt
2336  do i=is-nt,ie+1+nt
2337  fx(i,j) = dy(i,j)*sina_u(i,j)*(q(i-1,j,k)-q(i,j,k))*rdxc(i,j)
2338  enddo
2339  enddo
2340 
2341  if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 2)
2342  do j=js-nt,je+1+nt
2343  do i=is-nt,ie+nt
2344  fy(i,j) = dx(i,j)*sina_v(i,j)*(q(i,j-1,k)-q(i,j,k))*rdyc(i,j)
2345  enddo
2346  enddo
2347 
2348  if ( nt==0 ) then
2349  do j=js,je
2350  do i=is,ie
2351  qdt(i,j,k) = q(i,j,k) + damp*rarea(i,j)*(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))
2352  enddo
2353  enddo
2354  else
2355  do j=js-nt,je+nt
2356  do i=is-nt,ie+nt
2357  q(i,j,k) = q(i,j,k) + damp*rarea(i,j)*(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))
2358  enddo
2359  enddo
2360  endif
2361  enddo
2362 
2363  enddo
2364 
2365  end subroutine del2_scalar
2366 
2367 
2368 end module nwp_nudge_nlm_mod
fms_mod
Definition: fms.F90:20
nwp_nudge_nlm_mod::del_r
real(fvprc), parameter del_r
Definition: nwp_nudge_nlm.F90:109
time_manager_mod::get_date
subroutine, public get_date(time, year, month, day, hour, minute, second, tick, err_msg)
Definition: time_manager.F90:1666
fv_mp_nlm_mod
Definition: fv_mp_nlm_mod.F90:24
nwp_nudge_nlm_mod::bk0
real(fvprc), dimension(:), allocatable bk0
Definition: nwp_nudge_nlm.F90:33
constants_mod
Definition: constants.F90:24
constants_mod::radius
real, parameter, public radius
Radius of the Earth [m].
Definition: constants.F90:72
fv_arrays_nlm_mod::real8
integer, parameter real8
Definition: fv_arrays_nlm.F90:34
nwp_nudge_nlm_mod::nobs_tc
integer, dimension(max_storm) nobs_tc
Definition: nwp_nudge_nlm.F90:115
nwp_nudge_nlm_mod::ibtrack
logical ibtrack
Definition: nwp_nudge_nlm.F90:64
nwp_nudge_nlm_mod::jdc
integer, dimension(:,:), allocatable jdc
Definition: nwp_nudge_nlm.F90:50
nwp_nudge_nlm_mod::kstart
integer kstart
Definition: nwp_nudge_nlm.F90:92
nwp_nudge_nlm_mod::breed_slp
subroutine breed_slp(time, dt, npz, ak, bk, ps, phis, delp, ua, va, u_dt, v_dt, pt, q, nwat, zvir)
Definition: nwp_nudge_nlm.F90:1743
nwp_nudge_nlm_mod::ps_dat
real(fvprc), dimension(:,:,:), allocatable ps_dat
Definition: nwp_nudge_nlm.F90:51
nwp_nudge_nlm_mod::nf_t
integer nf_t
Definition: nwp_nudge_nlm.F90:89
nwp_nudge_nlm_mod::u_dat
real *4, dimension(:,:,:,:), allocatable u_dat
Definition: nwp_nudge_nlm.F90:52
nwp_nudge_nlm_mod::nobs_max
integer, parameter nobs_max
Definition: nwp_nudge_nlm.F90:112
nwp_nudge_nlm_mod::nstorms
integer nstorms
Definition: nwp_nudge_nlm.F90:114
nwp_nudge_nlm_mod::del2_uv
subroutine del2_uv(du, dv, cd, kmd, ntimes)
Definition: nwp_nudge_nlm.F90:2265
nwp_nudge_nlm_mod::get_obs
subroutine get_obs(Time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, u_obs, v_obs, t_obs, q_obs, tpw_dat, phis, gz_int, npz)
Definition: nwp_nudge_nlm.F90:460
nwp_nudge_nlm_mod::q_dat
real *4, dimension(:,:,:,:), allocatable q_dat
Definition: nwp_nudge_nlm.F90:52
nwp_nudge_nlm_mod::lon
real(fvprc), dimension(:), allocatable lon
Definition: nwp_nudge_nlm.F90:34
nwp_nudge_nlm_mod::get_int_hght
subroutine get_int_hght(h_int, npz, ak, bk, ps, delp, ps0, tv)
Definition: nwp_nudge_nlm.F90:1084
nwp_nudge_nlm_mod::nudge_debug
logical nudge_debug
Definition: nwp_nudge_nlm.F90:65
nwp_nudge_nlm_mod::nfile_total
integer nfile_total
Definition: nwp_nudge_nlm.F90:58
mpp_domains_mod::mpp_update_domains
Definition: mpp_domains.F90:1068
nwp_nudge_nlm_mod::ncep2fms
subroutine ncep2fms(sst)
Definition: nwp_nudge_nlm.F90:998
nwp_nudge_nlm_mod::slp_obs_init
subroutine slp_obs_init
Definition: nwp_nudge_nlm.F90:2042
nwp_nudge_nlm_mod::tau_t
real(fvprc) tau_t
Definition: nwp_nudge_nlm.F90:81
nwp_nudge_nlm_mod::nudge_winds
logical nudge_winds
Definition: nwp_nudge_nlm.F90:68
nwp_nudge_nlm_mod::get_ncep_analysis
subroutine get_ncep_analysis(ps, u, v, t, q, zvir, ts, nfile, fname)
Definition: nwp_nudge_nlm.F90:725
constants_mod::rdgas
real, parameter, public rdgas
Gas constant for dry air [J/kg/deg].
Definition: constants.F90:77
nwp_nudge_nlm_mod::ak0
real(fvprc), dimension(:), allocatable ak0
Definition: nwp_nudge_nlm.F90:33
nwp_nudge_nlm_mod::master
logical master
Definition: nwp_nudge_nlm.F90:37
nwp_nudge_nlm_mod::tau_hght
real(fvprc) tau_hght
Definition: nwp_nudge_nlm.F90:83
nwp_nudge_nlm_mod::x_obs
real *4, dimension(nobs_max, max_storm) x_obs
Definition: nwp_nudge_nlm.F90:116
fms_io_mod
Definition: fms_io.F90:20
nwp_nudge_nlm_mod::y_obs
real *4, dimension(nobs_max, max_storm) y_obs
Definition: nwp_nudge_nlm.F90:117
fv_arrays_nlm_mod::real4
integer, parameter real4
Definition: fv_arrays_nlm.F90:33
nwp_nudge_nlm_mod::nwp_nudge_init
subroutine, public nwp_nudge_init(npz, zvir, ak, bk, ts, phis)
Definition: nwp_nudge_nlm.F90:595
nwp_nudge_nlm_mod::mslp_obs
real *4, dimension(nobs_max, max_storm) mslp_obs
Definition: nwp_nudge_nlm.F90:118
nwp_nudge_nlm_mod::nfile
integer nfile
Definition: nwp_nudge_nlm.F90:43
nwp_nudge_nlm_mod::time_nudge
real(fvprc) time_nudge
Definition: nwp_nudge_nlm.F90:40
nwp_nudge_nlm_mod::leap_year
logical function leap_year(ny)
Definition: nwp_nudge_nlm.F90:2203
nwp_nudge_nlm_mod::version
character(len=128) version
Definition: nwp_nudge_nlm.F90:25
mpp_mod
Definition: mpp.F90:39
nwp_nudge_nlm_mod::track_file_name
character(len=128) track_file_name
Definition: nwp_nudge_nlm.F90:57
fms_mod::check_nml_error
integer function, public check_nml_error(IOSTAT, NML_NAME)
Definition: fms.F90:658
fv_grid_utils_nlm_mod
Definition: fv_grid_utils_nlm.F90:20
fv_timing_nlm_mod
Definition: fv_timing_nlm.F90:20
nwp_nudge_nlm_mod::time_tc
real *4, dimension(nobs_max, max_storm) time_tc
Definition: nwp_nudge_nlm.F90:121
nwp_nudge_nlm_mod::p_trop
real(fvprc) p_trop
Definition: nwp_nudge_nlm.F90:47
nwp_nudge_nlm_mod::time_varying
logical time_varying
Definition: nwp_nudge_nlm.F90:72
nwp_nudge_nlm_mod::max_storm
integer, parameter max_storm
Definition: nwp_nudge_nlm.F90:111
mpp_domains_mod
Definition: mpp_domains.F90:128
nwp_nudge_nlm_mod::kbot_q
integer kbot_q
Definition: nwp_nudge_nlm.F90:97
time_manager_mod
Definition: time_manager.F90:19
nwp_nudge_nlm_mod::r_inc
real(fvprc) r_inc
Definition: nwp_nudge_nlm.F90:108
fv_control_nlm_mod::npx
integer, pointer npx
Definition: fv_control_nlm.F90:141
nwp_nudge_nlm_mod::pmaxmin
subroutine pmaxmin(qname, a, imax, jmax, fac)
Definition: nwp_nudge_nlm.F90:2229
nwp_nudge_nlm_mod::tau_tpw
real(fvprc) tau_tpw
Definition: nwp_nudge_nlm.F90:79
nwp_nudge_nlm_mod::nudge_t
logical nudge_t
Definition: nwp_nudge_nlm.F90:66
field_exist
int field_exist(const char *file, const char *name)
Definition: read_mosaic.c:69
nwp_nudge_nlm_mod::module_is_initialized
logical module_is_initialized
Definition: nwp_nudge_nlm.F90:36
nwp_nudge_nlm_mod::id2
integer, dimension(:,:), allocatable id2
Definition: nwp_nudge_nlm.F90:50
nwp_nudge_nlm_mod::remap_tq
subroutine remap_tq(npz, ak, bk, ps, delp, t, q, kmd, ps0, ta, qa, zvir)
Definition: nwp_nudge_nlm.F90:1141
nwp_nudge_nlm_mod::do_nwp_nudge
subroutine, public do_nwp_nudge(Time, dt, npz, ps_dt, u_dt, v_dt, t_dt, q_dt, zvir, ak, bk, ts, ps, delp, ua, va, pt, nwat, q, phis)
Definition: nwp_nudge_nlm.F90:136
nwp_nudge_nlm_mod::v_dat
real *4, dimension(:,:,:,:), allocatable v_dat
Definition: nwp_nudge_nlm.F90:52
nwp_nudge_nlm_mod::tau_virt
real(fvprc) tau_virt
Definition: nwp_nudge_nlm.F90:82
nwp_nudge_nlm_mod::gz0
real(fvprc), dimension(:,:), allocatable gz0
Definition: nwp_nudge_nlm.F90:53
fms_io_mod::field_size
subroutine, public field_size(filename, fieldname, siz, field_found, domain, no_domain)
Definition: fms_io.F90:4941
fv_timing_nlm_mod::timing_on
subroutine timing_on(blk_name)
Definition: fv_timing_nlm.F90:101
nwp_nudge_nlm_mod::kord_data
integer kord_data
Definition: nwp_nudge_nlm.F90:60
nwp_nudge_nlm_mod::tc_mask
logical tc_mask
Definition: nwp_nudge_nlm.F90:62
nwp_nudge_nlm_mod::get_slp_obs
subroutine get_slp_obs(time, nobs, lon_obs, lat_obs, mslp, slp_out, r_out, time_obs, x_o, y_o, slp_o, r_vor, p_vor, stime, fact)
Definition: nwp_nudge_nlm.F90:1967
nwp_nudge_nlm_mod::k_breed
integer k_breed
Definition: nwp_nudge_nlm.F90:45
time_manager_mod::time_type
Definition: time_manager.F90:166
nwp_nudge_nlm_mod::kbot_winds
integer kbot_winds
Definition: nwp_nudge_nlm.F90:95
nwp_nudge_nlm_mod::del2_scalar
subroutine del2_scalar(qdt, cd, kmd, ntimes)
Definition: nwp_nudge_nlm.F90:2304
constants_mod::cp_air
real, parameter, public cp_air
Specific heat capacity of dry air at constant pressure [J/kg/deg].
Definition: constants.F90:83
nwp_nudge_nlm_mod::r_min
real(fvprc) r_min
Definition: nwp_nudge_nlm.F90:107
nwp_nudge_nlm_mod::s2c
real(fvprc), dimension(:,:,:), allocatable s2c
Definition: nwp_nudge_nlm.F90:49
fv_arrays_nlm_mod
Definition: fv_arrays_nlm.F90:20
nwp_nudge_nlm_mod::tau_q
real(fvprc) tau_q
Definition: nwp_nudge_nlm.F90:78
nwp_nudge_nlm_mod::breed_slp_inline
subroutine, public breed_slp_inline(nstep, dt, npz, ak, bk, phis, pe, pk, peln, delp, u, v, pt, q, nwat, zvir)
Definition: nwp_nudge_nlm.F90:1373
nwp_nudge_nlm_mod::nudge_q
logical nudge_q
Definition: nwp_nudge_nlm.F90:67
fv_grid_tools_nlm_mod
Definition: fv_grid_tools_nlm.F90:20
fv_grid_utils_nlm_mod::great_circle_dist
real function, public great_circle_dist(q1, q2, radius)
Definition: fv_grid_utils_nlm.F90:1968
nwp_nudge_nlm_mod::jm
integer jm
Definition: nwp_nudge_nlm.F90:31
nwp_nudge_nlm_mod::lat
real(fvprc), dimension(:), allocatable lat
Definition: nwp_nudge_nlm.F90:34
nwp_nudge_nlm_mod::k_trop
integer k_trop
Definition: nwp_nudge_nlm.F90:46
nwp_nudge_nlm_mod::q_rat
real(fvprc) q_rat
Definition: nwp_nudge_nlm.F90:86
nwp_nudge_nlm_mod::q_min
real(fvprc) q_min
Definition: nwp_nudge_nlm.F90:85
nwp_nudge_nlm_mod
Definition: nwp_nudge_nlm.F90:1
nwp_nudge_nlm_mod::deg2rad
real(fvprc) deg2rad
Definition: nwp_nudge_nlm.F90:39
constants_mod::grav
real, parameter, public grav
Acceleration due to gravity [m/s^2].
Definition: constants.F90:76
nwp_nudge_nlm_mod::nudge_virt
logical nudge_virt
Definition: nwp_nudge_nlm.F90:69
nwp_nudge_nlm_mod::no_obs
logical no_obs
Definition: nwp_nudge_nlm.F90:38
nwp_nudge_nlm_mod::p_wvp
real(fvprc) p_wvp
Definition: nwp_nudge_nlm.F90:59
nwp_nudge_nlm_mod::nf_uv
integer nf_uv
Definition: nwp_nudge_nlm.F90:88
fv_mapz_nlm_mod
Definition: fv_mapz_nlm.F90:23
fv_mp_nlm_mod::gid
integer, public gid
Definition: fv_mp_nlm_mod.F90:2464
nwp_nudge_nlm_mod::tau_winds
real(fvprc) tau_winds
Definition: nwp_nudge_nlm.F90:80
fv_mapz_nlm_mod::mappm
subroutine, public mappm(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord, ptop)
Definition: fv_mapz_nlm.F90:3225
mpp_mod::mpp_error
Definition: mpp.F90:402
nwp_nudge_nlm_mod::strong_mask
logical strong_mask
Definition: nwp_nudge_nlm.F90:63
fv_grid_utils_nlm_mod::intp_great_circle
subroutine, public intp_great_circle(beta, p1, p2, x_o, y_o)
Definition: fv_grid_utils_nlm.F90:1824
nwp_nudge_nlm_mod::year_track_data
integer year_track_data
Definition: nwp_nudge_nlm.F90:110
max
#define max(a, b)
Definition: mosaic_util.h:33
nwp_nudge_nlm_mod::time_track
logical time_track
Definition: nwp_nudge_nlm.F90:73
nwp_nudge_nlm_mod::breed_vortex
logical breed_vortex
Definition: nwp_nudge_nlm.F90:103
tp_core_nlm_mod
Definition: tp_core_nlm.F90:20
nwp_nudge_nlm_mod::nudge_tpw
logical nudge_tpw
Definition: nwp_nudge_nlm.F90:71
nwp_nudge_nlm_mod::rad2deg
real(fvprc) rad2deg
Definition: nwp_nudge_nlm.F90:39
nwp_nudge_nlm_mod::tagname
character(len=128) tagname
Definition: nwp_nudge_nlm.F90:26
fms_io_mod::read_data
Definition: fms_io.F90:323
tp_core_nlm_mod::copy_corners
subroutine, public copy_corners(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
Definition: tp_core_nlm.F90:216
nwp_nudge_nlm_mod::remap_coef
subroutine remap_coef
Definition: nwp_nudge_nlm.F90:925
fv_diagnostics_nlm_mod::prt_maxmin
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
Definition: fv_diagnostics_nlm.F90:2771
nwp_nudge_nlm_mod::get_tc_mask
subroutine get_tc_mask(time, mask)
Definition: nwp_nudge_nlm.F90:1329
nwp_nudge_nlm_mod::nudge_hght
logical nudge_hght
Definition: nwp_nudge_nlm.F90:70
nwp_nudge_nlm_mod::nfile_max
integer, parameter nfile_max
Definition: nwp_nudge_nlm.F90:42
nwp_nudge_nlm_mod::calday
real(fvprc) function calday(year, month, day, hour, minute, sec)
Definition: nwp_nudge_nlm.F90:2169
nwp_nudge_nlm_mod::slp_env
real(fvprc) slp_env
Definition: nwp_nudge_nlm.F90:106
fv_mp_nlm_mod::masterproc
integer, public masterproc
Definition: fv_mp_nlm_mod.F90:2463
nwp_nudge_nlm_mod::km
integer km
Definition: nwp_nudge_nlm.F90:32
nwp_nudge_nlm_mod::time_interval
integer time_interval
Definition: nwp_nudge_nlm.F90:41
nwp_nudge_nlm_mod::nwp_nudge_end
subroutine, public nwp_nudge_end
Definition: nwp_nudge_nlm.F90:1303
min
#define min(a, b)
Definition: mosaic_util.h:32
time_manager_mod::get_time
subroutine, public get_time(Time, seconds, days, ticks, err_msg)
Definition: time_manager.F90:493
fv_control_nlm_mod::npy
integer, pointer npy
Definition: fv_control_nlm.F90:142
fv_arrays_nlm_mod::fvprc
integer, parameter fvprc
Definition: fv_arrays_nlm.F90:38
nwp_nudge_nlm_mod::tau_vortex
real(fvprc) tau_vortex
Definition: nwp_nudge_nlm.F90:104
fv_diagnostics_nlm_mod
Definition: fv_diagnostics_nlm.F90:20
constants_mod::kappa
real, parameter, public kappa
RDGAS / CP_AIR [dimensionless].
Definition: constants.F90:82
nwp_nudge_nlm_mod::rad_out
real *4, dimension(nobs_max, max_storm) rad_out
Definition: nwp_nudge_nlm.F90:120
fv_diagnostics_nlm_mod::fv_time
type(time_type), public fv_time
Definition: fv_diagnostics_nlm.F90:61
nwp_nudge_nlm_mod::remap_uv
subroutine remap_uv(npz, ak, bk, ps, delp, u, v, kmd, ps0, u0, v0)
Definition: nwp_nudge_nlm.F90:1228
nwp_nudge_nlm_mod::kbot_t
integer kbot_t
Definition: nwp_nudge_nlm.F90:96
nwp_nudge_nlm_mod::id1
integer, dimension(:,:), allocatable id1
Definition: nwp_nudge_nlm.F90:50
fv_control_nlm_mod
Definition: fv_control_nlm.F90:25
nwp_nudge_nlm_mod::im
integer im
Definition: nwp_nudge_nlm.F90:30
crtm_parameters::pi
real(fp), parameter, public pi
Definition: CRTM_Parameters.f90:86
fv_timing_nlm_mod::timing_off
subroutine timing_off(blk_name)
Definition: fv_timing_nlm.F90:175
nwp_nudge_nlm_mod::file_names
character(len=128), dimension(nfile_max) file_names
Definition: nwp_nudge_nlm.F90:56
nwp_nudge_nlm_mod::t_dat
real *4, dimension(:,:,:,:), allocatable t_dat
Definition: nwp_nudge_nlm.F90:52
nwp_nudge_nlm_mod::mslp_out
real *4, dimension(nobs_max, max_storm) mslp_out
Definition: nwp_nudge_nlm.F90:119