FV3DYCORE  Version 1.1.0
fv_nudge.F90
Go to the documentation of this file.
1 
2 !***********************************************************************
3 !* GNU Lesser General Public License
4 !*
5 !* This file is part of the FV3 dynamical core.
6 !*
7 !* The FV3 dynamical core is free software: you can redistribute it
8 !* and/or modify it under the terms of the
9 !* GNU Lesser General Public License as published by the
10 !* Free Software Foundation, either version 3 of the License, or
11 !* (at your option) any later version.
12 !*
13 !* The FV3 dynamical core is distributed in the hope that it will be
14 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
15 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16 !* See the GNU General Public License for more details.
17 !*
18 !* You should have received a copy of the GNU Lesser General Public
19 !* License along with the FV3 dynamical core.
20 !* If not, see <http://www.gnu.org/licenses/>.
21 !***********************************************************************
22 #ifdef OVERLOAD_R4
23 #define _GET_VAR1 get_var1_real
24 #else
25 #define _GET_VAR1 get_var1_double
26 #endif
27 
31 
33 
34 ! <table>
35 ! <tr>
36 ! <th>Module Name</th>
37 ! <th>Functions Included</th>
38 ! </tr>
39 ! <tr>
40 ! <td>constants_mod</td>
41 ! <td>pi=>pi_8, grav, rdgas, cp_air, kappa, cnst_radius =>radius</td>
42 ! </tr>
43 ! <tr>
44 ! <td>external_sst_mod</td>
45 ! <td>i_sst, j_sst, sst_ncep, sst_anom, forecast_mode</td>
46 ! </tr>
47 ! <tr>
48 ! <td>diag_manager_mod</td>
49 ! <td>register_diag_field, send_data</td>
50 ! </tr>
51 ! <tr>
52 ! <td>fms_mod</td>
53 ! <td>write_version_number, open_namelist_file, check_nml_error,
54 ! file_exist, close_file</td>
55 ! </tr>
56 ! <tr>
57 ! <td>fv_arrays_mod</td>
58 ! <td>fv_grid_type, fv_grid_bounds_type, fv_nest_type, R_GRID</td>
59 ! </tr>
60 ! <tr>
61 ! <td>fv_diagnostics_mod</td>
62 ! <td>prt_maxmin, fv_time</td>
63 ! </tr>
64 ! <tr>
65 ! <td>fv_grid_utils_mod</td>
66 ! <td>great_circle_dist, intp_great_circle,
67 ! latlon2xyz, vect_cross, normalize_vect</td>
68 ! </tr>
69 ! <tr>
70 ! <td>fv_mp_mod</td>
71 ! <td>nmp_reduce_sum, mp_reduce_min, mp_reduce_max, is_master</td>
72 ! </tr>
73 ! <tr>
74 ! <td>fv_mapz_mod</td>
75 ! <td>mappm</td>
76 ! </tr>
77 ! <tr>
78 ! <td>fv_timing_mod</td>
79 ! <td>timing_on, timing_off</td>
80 ! </tr>
81 ! <tr>
82 ! <td>mpp_mod</td>
83 ! <td>mpp_error, FATAL, stdlog, get_unit, mpp_pe</td>
84 ! </tr>
85 ! <tr>
86 ! <td>sim_nc_mod</td>
87 ! <td>open_ncfile, close_ncfile, get_ncdim1, get_var1_double,
88 ! get_var3_r4, get_var1_real</td>
89 ! </tr>
90 ! <tr>
91 ! <td>mpp_domains_mod</td>
92 ! <td>mpp_update_domains, domain2d</td>
93 ! </tr>
94 ! <tr>
95 ! <td>time_manager_mod</td>
96 ! <td>time_type, get_time, get_date</td>
97 ! </tr>
98 ! <tr>
99 ! <td>tp_core_mod</td>
100 ! <td>copy_corners</td>
101 ! </tr>
102 ! </table>
103 
104  use external_sst_mod, only: i_sst, j_sst, sst_ncep, sst_anom, forecast_mode
105  use diag_manager_mod, only: register_diag_field, send_data
106  use constants_mod, only: pi=>pi_8, grav, rdgas, cp_air, kappa, cnst_radius =>radius
107  use fms_mod, only: write_version_number, open_namelist_file, &
108  check_nml_error, file_exist, close_file
109 !use fms_io_mod, only: field_size
110  use mpp_mod, only: mpp_error, fatal, stdlog, get_unit, mpp_pe
111  use mpp_domains_mod, only: mpp_update_domains, domain2d
112  use time_manager_mod, only: time_type, get_time, get_date
113 
117  use tp_core_mod, only: copy_corners
118  use fv_mapz_mod, only: mappm
119  use fv_mp_mod, only: mp_reduce_sum, mp_reduce_min, mp_reduce_max, is_master
121 
125 #ifdef MULTI_GASES
126  use multi_gases_mod, only: virq, virqd, vicpqd, num_gas
127 #endif
128 
129  implicit none
130  private
131 
132  real(kind=R_GRID), parameter :: radius = cnst_radius
133 
134 ! version number of this module
135 ! Include variable "version" to be written to log file.
136 #include<file_version.h>
137 
138  logical :: do_adiabatic_init
139 
141  public do_adiabatic_init
142  integer im
143  integer jm
144  integer km
145  real, allocatable:: ak0(:), bk0(:)
146  real, allocatable:: lat(:), lon(:)
147 
148  logical :: module_is_initialized = .false.
149  logical :: master
150  logical :: no_obs
151  real :: deg2rad, rad2deg
152  real :: time_nudge = 0.
153  integer :: time_interval = 6*3600
154 ! ---> h1g, enhance the max. analysis data files, 2012-10-22
155 ! integer, parameter :: nfile_max = 125
156  integer, parameter :: nfile_max = 29280
157 ! <--- h1g, 2012-10-22
158  integer :: nfile
159 
160  integer :: k_breed = 0
161  integer :: k_trop = 0
162  real :: p_trop = 950.e2
163  real :: dps_min = 50.
164  real :: del2_cd = 0.16
165 
166  real, allocatable:: s2c(:,:,:)
167  integer, allocatable:: id1(:,:), id2(:,:), jdc(:,:)
168  real, allocatable :: ps_dat(:,:,:)
169  real(KIND=4), allocatable, dimension(:,:,:,:):: u_dat, v_dat, t_dat, q_dat
170  real(KIND=4), allocatable, dimension(:,:,:):: gz3
171  real, allocatable:: gz0(:,:)
172 
173 ! Namelist variables:
174 ! ---> h1g, add the list of input NCEP analysis data files, 2012-10-22
175  character(len=128):: input_fname_list =""
176  character(len=128):: analysis_file_first =""
178  character(len=128):: analysis_file_last=""
180 
181  real :: p_relax = 30.e2
183 
184  real :: p_norelax = 0.0
185 ! <--- h1g, 2012-10-22
186 
187  character(len=128):: file_names(nfile_max)
188  character(len=128):: track_file_name
189  integer :: nfile_total = 0
190  real :: p_wvp = 100.e2
191  integer :: kord_data = 8
192 
193  real :: mask_fac = 0.25
194 
195  logical :: t_is_tv = .false.
196  logical :: use_pt_inc = .false.
197  logical :: use_high_top = .false.
198  logical :: add_bg_wind = .true.
199  logical :: conserve_mom = .true.
200  logical :: conserve_hgt = .true.
201  logical :: tc_mask = .false.
202  logical :: strong_mask = .false.
203  logical :: ibtrack = .true.
204  logical :: nudge_debug = .false.
205  logical :: do_ps_bias = .false.
206  logical :: nudge_ps = .false.
207  logical :: nudge_q = .false.
208  logical :: nudge_winds = .true.
209  logical :: nudge_virt = .true.
210  logical :: nudge_hght = .true.
211  logical :: time_varying = .true.
212  logical :: print_end_breed = .true.
213  logical :: print_end_nudge = .true.
214 
215 !Nudging time-scales (seconds):
216 !note, however, the effective time-scale is 2X smaller (stronger) due
217 ! to the use of the time-varying weighting factor
218  real :: tau_ps = 21600.
219  real :: tau_q = 86400.
220  real :: tau_winds = 21600.
221  real :: tau_virt = 43200.
222  real :: tau_hght = 43200.
223 
224  real :: q_min = 1.e-8
225 
226  integer :: jbeg, jend
227  integer :: nf_uv = 0
228  integer :: nf_ps = 0
229  integer :: nf_t = 2
230  integer :: nf_ht = 1
231 
232 ! starting layer (top layer is sponge layer and is skipped)
233  integer :: kstart = 2
234 
235 ! skip "kbot" layers
236  integer :: kbot_winds = 0
237  integer :: kbot_t = 0
238  integer :: kbot_q = 0
239  logical :: analysis_time
240 
241 !-- Tropical cyclones --------------------------------------------------------------------
242 
243 ! track dataset: 'INPUT/tropical_cyclones.txt'
244 
245  logical :: breed_srf_w = .false.
246  real :: grid_size = 28.e3
247  real :: tau_vt_slp = 1200.
248  real :: tau_vt_wind = 1200.
249  real :: tau_vt_rad = 4.0
250 
251  real :: pt_lim = 0.2
252  real :: slp_env = 101010.
253  real :: pre0_env = 100000.
254  real, parameter:: tm_max = 315.
255 !------------------
256  real:: r_lo = 2.0
257  real:: r_hi = 5.0 ! try 4.0?
258 !------------------
259  real:: r_fac = 1.2
260  real :: r_min = 200.e3
261  real :: r_inc = 25.e3
262  real, parameter:: del_r = 50.e3
263  real:: elapsed_time = 0.0
264  real:: nudged_time = 1.e12
265  ! usage example: set to 43200. to do inline vortex breeding
266  ! for only the first 12 hours
267  ! In addition, specify only 3 analysis files (12 hours)
268  integer:: year_track_data
269  integer, parameter:: max_storm = 140
270  integer, parameter:: nobs_max = 125
271 
272  integer :: nstorms = 0
273  integer :: nobs_tc(max_storm)
274  integer :: min_nobs = 16
275  real :: min_mslp = 1009.e2
276  real(KIND=4):: x_obs(nobs_max,max_storm)
277  real(KIND=4):: y_obs(nobs_max,max_storm)
278  real(KIND=4):: wind_obs(nobs_max,max_storm)
279  real(KIND=4):: mslp_obs(nobs_max,max_storm)
280  real(KIND=4):: mslp_out(nobs_max,max_storm)
281  real(KIND=4):: rad_out(nobs_max,max_storm)
282  real(KIND=4):: time_tc(nobs_max,max_storm)
283 !------------------------------------------------------------------------------------------
284  integer :: id_ht_err
285 
286  !!! CLEANUP Module pointers
287  integer :: is, ie, js, je
288  integer :: isd, ied, jsd, jed
289 
290  namelist /fv_nwp_nudge_nml/t_is_tv, nudge_ps, nudge_virt, nudge_hght, nudge_q, nudge_winds, &
297  input_fname_list, analysis_file_first, analysis_file_last, p_relax, p_norelax !h1g, add 3 namelist variables, 2012-20-22
298 
299  contains
300 
303  subroutine fv_nwp_nudge ( Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt, zvir, ptop, &
304  ak, bk, ts, ps, delp, ua, va, pt, nwat, q, phis, gridstruct, &
305  bd, domain )
307  type(time_type), intent(in):: Time
308  integer, intent(IN):: npx, npy
309  integer, intent(in):: npz
310  integer, intent(in):: nwat
311  real, intent(in):: dt
312  real, intent(in):: zvir, ptop
313  type(domain2d), intent(INOUT), target :: domain
314  type(fv_grid_bounds_type), intent(IN) :: bd
315  real, intent(in ), dimension(npz+1):: ak, bk
316  real, intent(in ), dimension(isd:ied,jsd:jed ):: phis
317  real, intent(inout), dimension(isd:ied,jsd:jed,npz):: pt, ua, va, delp
318 ! pt as input is true tempeature
319 #ifdef MULTI_GASES
320  real, intent(inout):: q(isd:ied,jsd:jed,npz,*)
321 #else
322  real, intent(inout):: q(isd:ied,jsd:jed,npz,nwat)
323 #endif
324  real, intent(inout), dimension(isd:ied,jsd:jed):: ps
325 ! Accumulated tendencies
326  real, intent(inout), dimension(isd:ied,jsd:jed,npz):: u_dt, v_dt
327  real, intent(out), dimension(is:ie,js:je,npz):: t_dt, q_dt
328  real, intent(out), dimension(is:ie,js:je):: ps_dt, ts
329 
330  type(fv_grid_type), intent(INOUT), target :: gridstruct
331 ! local:
332  real:: h2(is:ie,js:je)
333  real:: m_err(is:ie,js:je)
334  real:: slp_n(is:ie,js:je)
335  real:: slp_m(is:ie,js:je)
336  real:: mask(is:ie,js:je)
337  real:: tv(is:ie,js:je)
338  real:: peln(is:ie,npz+1)
339  real:: pe2(is:ie, npz+1)
340  real:: ptmp
341  real, allocatable :: ps_obs(:,:)
342  real, allocatable, dimension(:,:,:):: u_obs, v_obs, t_obs, q_obs
343  real, allocatable, dimension(:,:,:):: du_obs, dv_obs
344  real:: ps_fac(is:ie,js:je)
345  integer :: seconds, days
346  integer :: i,j,k, iq, kht, n
347  real :: factor, rms, bias, co
348  real :: rdt, press(npz), profile(npz), prof_t(npz), prof_q(npz), du, dv
349  logical used
350 
351 
352  real(kind=R_GRID), pointer, dimension(:,:,:) :: agrid
353  real, pointer, dimension(:,:) :: rarea, area
354 
355  real, pointer, dimension(:,:) :: sina_u, sina_v
356  real, pointer, dimension(:,:,:) :: sin_sg
357  real(kind=R_GRID), pointer, dimension(:,:,:) :: vlon, vlat
358 
359  real, pointer, dimension(:,:) :: dx, dy, rdxc, rdyc
360 
361  real(kind=R_GRID), pointer :: da_min
362 
363  logical, pointer :: nested, sw_corner, se_corner, nw_corner, ne_corner
364 
365  if ( .not. module_is_initialized ) then
366  call mpp_error(fatal,'==> Error from fv_nwp_nudge: module not initialized')
367  endif
368  agrid => gridstruct%agrid_64
369  area => gridstruct%area
370  rarea => gridstruct%rarea
371 
372  vlon => gridstruct%vlon
373  vlat => gridstruct%vlat
374  sina_u => gridstruct%sina_u
375  sina_v => gridstruct%sina_v
376  sin_sg => gridstruct%sin_sg
377 
378  dx => gridstruct%dx
379  dy => gridstruct%dy
380  rdxc => gridstruct%rdxc
381  rdyc => gridstruct%rdyc
382 
383  da_min => gridstruct%da_min
384 
385  nested => gridstruct%nested
386  sw_corner => gridstruct%sw_corner
387  se_corner => gridstruct%se_corner
388  nw_corner => gridstruct%nw_corner
389  ne_corner => gridstruct%ne_corner
390 
391  if ( no_obs ) then
392 #ifndef DYCORE_SOLO
393  forecast_mode = .true.
394 #endif
395  return
396  endif
397 
398  call get_time (time, seconds, days)
399 
400  do j=js,je
401  do i=is,ie
402  mask(i,j) = 1.
403  enddo
404  enddo
405  if ( tc_mask ) call get_tc_mask(time, mask, agrid)
406 
407 ! The following profile is suitable only for nwp purposes; if the analysis has a good representation
408 ! of the strat-meso-sphere the profile for upper layers should be changed.
409 
410  profile(:) = 1.
411 
412 !$OMP parallel do default(none) shared(npz,press,ak,bk,P_relax,P_norelax,profile)
413  do k=1,npz
414  press(k) = 0.5*(ak(k) + ak(k+1)) + 0.5*(bk(k)+bk(k+1))*1.e5
415  if ( press(k) < p_relax ) then
416  profile(k) = max(0.01, press(k)/p_relax)
417  endif
418 
419  ! above P_norelax, no nudging. added by h1g
420  if( press(k) < p_norelax ) profile(k) = 0.0
421  enddo
422  profile(1) = 0.
423 
424 ! Thermodynamics:
425  prof_t(:) = 1.
426 !$OMP parallel do default(none) shared(npz,press,prof_t)
427  do k=1,npz
428  if ( press(k) < 10.e2 ) then
429  prof_t(k) = max(0.01, press(k)/10.e2)
430  endif
431  enddo
432  prof_t(1) = 0.
433 
434 ! Water vapor:
435  prof_q(:) = 1.
436 !$OMP parallel do default(none) shared(npz,press,prof_q)
437  do k=1,npz
438  if ( press(k) < 300.e2 ) then
439  prof_q(k) = max(0., press(k)/300.e2)
440  endif
441  enddo
442  prof_q(1) = 0.
443 
444 ! Height
445  if ( k_trop == 0 ) then
446  k_trop = 2
447  do k=2,npz-1
448  ptmp = ak(k+1) + bk(k+1)*1.e5
449  if ( ptmp > p_trop ) then
450  k_trop = k
451  exit
452  endif
453  enddo
454  endif
455 
456  if ( nudge_virt .and. nudge_hght ) then
457  kht = k_trop
458  else
459  kht = npz-kbot_t
460  endif
461 
462  if ( time_varying ) then
463  factor = 1. + cos(real(mod(seconds,time_interval))/real(time_interval)*2.*pi)
464  factor = max(1.e-5, factor)
465  else
466  factor = 1.
467  endif
468 
469  if ( do_adiabatic_init ) factor = 2.*factor
470 
471  allocate (ps_obs(is:ie,js:je) )
472  allocate ( t_obs(is:ie,js:je,npz) )
473  allocate ( q_obs(is:ie,js:je,npz) )
474 
475  if ( nudge_winds ) then
476  allocate (u_obs(is:ie,js:je,npz) )
477  allocate (v_obs(is:ie,js:je,npz) )
478  endif
479 
480 
481  call get_obs(time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, pt, nwat, q, u_obs, v_obs, t_obs, q_obs, &
482  phis, ua, va, u_dt, v_dt, npx, npy, npz, factor, mask, ptop, bd, gridstruct, domain)
483 ! *t_obs* is virtual temperature
484 
485  if ( no_obs ) then
486  deallocate (ps_obs)
487  deallocate (t_obs)
488  deallocate (q_obs)
489  if ( nudge_winds ) then
490  deallocate (u_obs)
491  deallocate (v_obs)
492  endif
493 #ifndef DYCORE_SOLO
494  forecast_mode = .true.
495 #endif
496  return
497  endif
498 
499  do j=js,je
500  do i=is,ie
501  if ( abs(ps(i,j)-ps_obs(i,j)) > 2.e2 ) then
502  ps_fac(i,j) = 2.e2 / abs(ps(i,j)-ps_obs(i,j))
503  else
504  ps_fac(i,j) = 1.
505  endif
506  enddo
507  enddo
508 
509  if( analysis_time ) then
510 !-------------------------------------------
511 ! Compute RMSE, bias, and correlation of SLP
512 !-------------------------------------------
513  do j=js,je
514  do i=is,ie
515 #ifdef MULTI_GASES
516  tv(i,j) = pt(i,j,npz)*virq(q(i,j,npz,1:num_gas))
517 #else
518  tv(i,j) = pt(i,j,npz)*(1.+zvir*q(i,j,npz,1))
519 #endif
520  enddo
521  enddo
522  call compute_slp(is, ie, js, je, tv, ps(is:ie,js:je), phis(is:ie,js:je), slp_m)
523  call compute_slp(is, ie, js, je, t_obs(is,js,npz), ps_obs, gz0, slp_n)
524 
525  if ( nudge_debug ) then
526  if(master) write(*,*) 'kht=', kht
527  call prt_maxmin('PS_o', ps_obs, is, ie, js, je, 0, 1, 0.01)
528  ptmp = 0.01*g0_sum(ps_obs, is, ie, js, je, 1, .false., isd, ied, jsd, jed, area)
529  if(master) write(*,*) 'Mean PS_o=', ptmp
530  call prt_maxmin('SLP_m', slp_m, is, ie, js, je, 0, 1, 0.01)
531  call prt_maxmin('SLP_o', slp_n, is, ie, js, je, 0, 1, 0.01)
532  endif
533 
534 !$OMP parallel do default(none) shared(is,ie,js,je,phis,gz0,m_err,mask,slp_m,slp_n)
535  do j=js,je
536  do i=is,ie
537  if ( abs(phis(i,j)-gz0(i,j))/grav > 2. ) then
538  m_err(i,j) = mask(i,j)*(slp_m(i,j) - slp_n(i,j))*2.*grav/abs(phis(i,j)-gz0(i,j))
539  else
540  m_err(i,j) = mask(i,j)*(slp_m(i,j) - slp_n(i,j))
541  endif
542  enddo
543  enddo
544 
545  call rmse_bias(m_err, rms, bias, area)
546  call corr(slp_m, slp_n, co, area)
547 
548  call prt_maxmin('SLP Error (mb)=', m_err, is, ie, js, je, 0, 1, 0.01)
549  if(master) write(*,*) 'SLP (Pa): RMS=', rms, ' Bias=', bias
550  if(master) write(*,*) 'SLP correlation=',co
551  endif
552 
553  if ( nudge_winds ) then
554 
555  allocate (du_obs(is:ie,js:je,npz) )
556  allocate (dv_obs(is:ie,js:je,npz) )
557 
558  du_obs = 0.
559  dv_obs = 0.
560 
561 ! Compute tendencies:
562  rdt = 1. / (tau_winds/factor + dt)
563 !$OMP parallel do default(none) shared(kstart,npz,kbot_winds,is,ie,js,je,du_obs,dv_obs, &
564 !$OMP profile,u_obs,v_obs,ua,va,rdt)
565  do k=kstart, npz - kbot_winds
566  do j=js,je
567  do i=is,ie
568  du_obs(i,j,k) = profile(k)*(u_obs(i,j,k)-ua(i,j,k))*rdt
569  dv_obs(i,j,k) = profile(k)*(v_obs(i,j,k)-va(i,j,k))*rdt
570  enddo
571  enddo
572  enddo
573 
574  if ( nf_uv>0 ) call del2_uv(du_obs, dv_obs, del2_cd, npz, nf_uv, bd, npx, npy, gridstruct, domain)
575 
576 !$OMP parallel do default(none) shared(kstart,kbot_winds,npz,is,ie,js,je,du_obs,dv_obs, &
577 !$OMP mask,ps_fac,u_dt,v_dt,ua,va,dt)
578  do k=kstart, npz - kbot_winds
579  if ( k==npz ) then
580  do j=js,je
581  do i=is,ie
582  du_obs(i,j,k) = du_obs(i,j,k) * mask(i,j) * ps_fac(i,j)
583  dv_obs(i,j,k) = dv_obs(i,j,k) * mask(i,j) * ps_fac(i,j)
584 !
585  u_dt(i,j,k) = u_dt(i,j,k) + du_obs(i,j,k)
586  v_dt(i,j,k) = v_dt(i,j,k) + dv_obs(i,j,k)
587  ua(i,j,k) = ua(i,j,k) + du_obs(i,j,k)*dt
588  va(i,j,k) = va(i,j,k) + dv_obs(i,j,k)*dt
589  enddo
590  enddo
591  else
592  do j=js,je
593  do i=is,ie
594 ! Apply TC mask
595  du_obs(i,j,k) = du_obs(i,j,k) * mask(i,j)
596  dv_obs(i,j,k) = dv_obs(i,j,k) * mask(i,j)
597 !
598  u_dt(i,j,k) = u_dt(i,j,k) + du_obs(i,j,k)
599  v_dt(i,j,k) = v_dt(i,j,k) + dv_obs(i,j,k)
600  ua(i,j,k) = ua(i,j,k) + du_obs(i,j,k)*dt
601  va(i,j,k) = va(i,j,k) + dv_obs(i,j,k)*dt
602  enddo
603  enddo
604  endif
605  enddo
606  endif
607 
608 !$OMP parallel do default(none) shared(is,ie,js,je,npz,t_dt)
609  do k=1,npz
610  do j=js,je
611  do i=is,ie
612  t_dt(i,j,k) = 0.
613  enddo
614  enddo
615  enddo
616 
617  if ( nudge_virt ) then
618  rdt = 1./(tau_virt/factor + dt)
619 !$OMP parallel do default(none) shared(is,ie,js,je,npz,kstart,kht,t_dt,prof_t,t_obs,zvir, &
620 #ifdef MULTI_GASES
621 !$OMP num_gas, &
622 #endif
623 !$OMP q,pt,rdt,ps_fac)
624  do k=kstart, kht
625  if ( k==npz ) then
626  do j=js,je
627  do i=is,ie
628 #ifdef MULTI_GASES
629  t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)/virq(q(i,j,k,1:num_gas))-pt(i,j,k))*rdt*ps_fac(i,j)
630 #else
631  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*ps_fac(i,j)
632 #endif
633  enddo
634  enddo
635  else
636  do j=js,je
637  do i=is,ie
638 #ifdef MULTI_GASES
639  t_dt(i,j,k) = prof_t(k)*(t_obs(i,j,k)/virq(q(i,j,k,1:num_gas))-pt(i,j,k))*rdt
640 #else
641  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
642 #endif
643  enddo
644  enddo
645  endif
646  enddo
647  endif
648 
649  if ( nudge_hght .and. kht<npz ) then ! averaged (in log-p) temperature
650  rdt = 1. / (tau_hght/factor + dt)
651 !$OMP parallel do default(none) shared(is,ie,js,je,npz,ak,h2,delp,kht,t_obs,pt,zvir, &
652 #ifdef MULTI_GASES
653 !$OMP num_gas, &
654 #endif
655 !$OMP q,rdt,ps_fac,mask,t_dt) &
656 !$OMP private(pe2, peln )
657  do j=js,je
658  do i=is,ie
659  pe2(i,1) = ak(1)
660  h2(i,j) = 0.
661  enddo
662  do k=1, npz
663  do i=is,ie
664  pe2(i,k+1) = pe2(i,k) + delp(i,j,k)
665  enddo
666  enddo
667  do k=kht+1, npz+1
668  do i=is,ie
669  peln(i,k) = log(pe2(i,k))
670  enddo
671  enddo
672 
673  do k=kht+1, npz
674  do i=is,ie
675 ! Difference between "Mean virtual tempearture" (pseudo height)
676 #ifdef MULTI_GASES
677  h2(i,j) = h2(i,j) + (t_obs(i,j,k)-pt(i,j,k)*virq(q(i,j,k,1:num_gas)))*(peln(i,k+1)-peln(i,k))
678 #else
679  h2(i,j) = h2(i,j) + (t_obs(i,j,k)-pt(i,j,k)*(1.+zvir*q(i,j,k,1)))*(peln(i,k+1)-peln(i,k))
680 #endif
681  enddo
682  enddo
683  do i=is,ie
684  h2(i,j) = h2(i,j) / (peln(i,npz+1)-peln(i,kht+1))
685  h2(i,j) = rdt*ps_fac(i,j)*h2(i,j)*mask(i,j)
686  enddo
687 
688  do k=kht+1, npz
689  do i=is,ie
690 #ifdef MULTI_GASES
691  t_dt(i,j,k) = h2(i,j) / virq(q(i,j,k,1:num_gas))
692 #else
693  t_dt(i,j,k) = h2(i,j) / (1.+zvir*q(i,j,k,1))
694 #endif
695  enddo
696  enddo
697  enddo ! j-loop
698  if(nudge_debug) call prt_maxmin('H2 increment=', h2, is, ie, js, je, 0, 1, dt)
699  endif
700 
701  if ( nudge_virt .or. nudge_hght ) then
702  if ( nf_t>0 ) call del2_scalar(t_dt, del2_cd, npz, nf_t, bd, npx, npy, gridstruct, domain)
703 !$OMP parallel do default(none) shared(kstart,npz,is,ie,js,je,pt,t_dt,dt,mask)
704  do k=kstart, npz
705  do j=js,je
706  do i=is,ie
707  pt(i,j,k) = pt(i,j,k) + t_dt(i,j,k)*dt*mask(i,j)
708  enddo
709  enddo
710  enddo
711  endif
712 
713  q_dt(:,:,:) = 0.
714  if ( nudge_q ) then
715  rdt = 1./(tau_q/factor + dt)
716 !$OMP parallel do default(none) shared(kstart,npz,kbot_q,is,ie,js,je,press,p_wvp,nwat, &
717 !$OMP q,delp,q_dt,prof_q,q_min,q_obs,rdt,mask,dt)
718  do k=kstart, npz - kbot_q
719  if ( press(k) > p_wvp ) then
720  do iq=2,nwat
721  do j=js,je
722  do i=is,ie
723  q(i,j,k,iq) = q(i,j,k,iq)*delp(i,j,k)
724  enddo
725  enddo
726  enddo
727 ! Specific humidity:
728  do j=js,je
729  do i=is,ie
730  delp(i,j,k) = delp(i,j,k)*(1.-q(i,j,k,1))
731  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)
732  q(i,j,k,1) = q(i,j,k,1) + q_dt(i,j,k)*dt
733  delp(i,j,k) = delp(i,j,k)/(1.-q(i,j,k,1))
734  enddo
735  enddo
736  do iq=2,nwat
737  do j=js,je
738  do i=is,ie
739  q(i,j,k,iq) = q(i,j,k,iq)/delp(i,j,k)
740  enddo
741  enddo
742  enddo
743  endif
744  enddo
745  endif
746 
747  ps_dt(:,:) = 0.
748 
749  deallocate ( t_obs )
750  deallocate ( q_obs )
751  deallocate ( ps_obs )
752 
753  if ( breed_srf_w .and. nudge_winds ) &
754  call breed_srf_winds(time, dt, npz, u_obs, v_obs, ak, bk, ps, phis, delp, ua, va, u_dt, v_dt, pt, q, nwat, zvir, gridstruct)
755 
756  if ( nudge_debug) then
757  call prt_maxmin('T increment=', t_dt, is, ie, js, je, 0, npz, dt)
758  call prt_maxmin('U increment=', du_obs, is, ie, js, je, 0, npz, dt)
759  call prt_maxmin('V increment=', dv_obs, is, ie, js, je, 0, npz, dt)
760  endif
761 
762  if ( nudge_winds ) then
763  deallocate ( u_obs )
764  deallocate ( v_obs )
765  deallocate (du_obs)
766  deallocate (dv_obs)
767  endif
768 
769  nullify(agrid)
770  nullify(area)
771  nullify(rarea)
772 
773  nullify(vlon)
774  nullify(vlat)
775  nullify(sina_u)
776  nullify(sina_v)
777  nullify(sin_sg)
778 
779  nullify(dx)
780  nullify(dy)
781  nullify(rdxc)
782  nullify(rdyc)
783 
784  nullify(da_min)
785 
786  nullify(nested)
787  nullify(sw_corner)
788  nullify(se_corner)
789  nullify(nw_corner)
790  nullify(ne_corner)
791 
792  end subroutine fv_nwp_nudge
793 
794 
795  subroutine ps_nudging(dt, factor, npz, ak, bk, ps_obs, mask, tm, ps, phis, delp, ua, va, pt, nwat, q, bd, npx, npy, gridstruct, domain)
796 ! Input
797  real, intent(in):: dt, factor
798  integer, intent(in):: npz, nwat, npx, npy
799  real, intent(in), dimension(npz+1):: ak, bk
800  type(fv_grid_bounds_type), intent(IN) :: bd
801  real, intent(in):: phis(isd:ied,jsd:jed)
802  real, intent(in), dimension(is:ie,js:je):: ps_obs, mask, tm
803  type(fv_grid_type), intent(IN), target :: gridstruct
804  type(domain2d), intent(INOUT) :: domain
805 ! Input/Output
806  real, intent(inout), dimension(isd:ied,jsd:jed):: ps
807  real, intent(inout), dimension(isd:ied,jsd:jed,npz):: delp, pt, ua, va
808  real, intent(inout):: q(isd:ied,jsd:jed,npz,nwat)
809 ! local
810  real, dimension(is:ie,js:je):: ps_dt
811  integer, parameter:: kmax = 100
812  real:: pn0(kmax+1), pk0(kmax+1), pe0(kmax+1)
813  real, dimension(is:ie,npz+1):: pe2, peln
814  real:: pst, dbk, pt0, rdt, bias
815  integer i, j, k, iq
816 
817  real, pointer, dimension(:,:) :: area
818 #ifdef MULTI_GASES
819  real :: kappax(km)
820  real :: pkx
821 #endif
822 
823  area => gridstruct%area
824 
825 ! Adjust interpolated ps to model terrain
826  if ( kmax < km ) call mpp_error(fatal,'==> KMAX must be larger than km')
827 
828  do j=js,je
829  do 666 i=is,ie
830 #ifdef MULTI_GASES
831  do k=1,km
832  kappax(k)= virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))
833  enddo
834 #endif
835  do k=1, km+1
836  pe0(k) = ak0(k) + bk0(k)*ps_obs(i,j)
837  pk0(k) = (ak0(k) + bk0(k)*ps_obs(i,j))**kappa
838  enddo
839  if( phis(i,j)>gz0(i,j) ) then
840  do k=km,1,-1
841  if( phis(i,j)<gz3(i,j,k) .and. phis(i,j) >= gz3(i,j,k+1) ) then
842  pst = pk0(k) + (pk0(k+1)-pk0(k))*(gz3(i,j,k)-phis(i,j))/(gz3(i,j,k)-gz3(i,j,k+1))
843  go to 666
844  endif
845  enddo
846  else
847  pn0(km ) = log(ak0(km) + bk0(km)*ps_obs(i,j))
848  pn0(km+1) = log(ps_obs(i,j))
849 ! Extrapolation into the ground using only the lowest layer potential temp
850 #ifdef MULTI_GASES
851  pkx = (pk0(km+1)-pk0(km))/(kappa*(pn0(km+1)-pn0(km)))
852  pt0 = tm(i,j)/exp(kappax(km)*log(pkx))
853  pkx = exp((kappax(km)-1.0)*log(pkx))
854  pst = pk0(km+1) + (gz0(i,j)-phis(i,j))/(cp_air*pt0*pkx)
855 #else
856  pt0 = tm(i,j)/(pk0(km+1)-pk0(km))*(kappa*(pn0(km+1)-pn0(km)))
857  pst = pk0(km+1) + (gz0(i,j)-phis(i,j))/(cp_air*pt0)
858 #endif
859  endif
860 #ifdef MULTI_GASES
861  k=km
862 666 ps_dt(i,j) = pst**(1./(kappa*kappax(k))) - ps(i,j)
863 #else
864 666 ps_dt(i,j) = pst**(1./kappa) - ps(i,j)
865 #endif
866  enddo ! j-loop
867 
868  if( nf_ps>0 ) call del2_scalar(ps_dt, del2_cd, 1, nf_ps, bd, npx, npy, gridstruct, domain)
869 
870  do j=js,je
871  do i=is,ie
872 ! Cap the errors:
873  ps_dt(i,j) = sign( min(10.e2, abs(ps_dt(i,j))), ps_dt(i,j) )
874  ps_dt(i,j) = mask(i,j)*ps_dt(i,j) * max( 0., 1.-abs(gz0(i,j)-phis(i,j))/(grav*500.) )
875  enddo
876  enddo
877 
878  if( do_ps_bias ) call ps_bias_correction( ps_dt, is, ie, js, je, isd, ied, jsd, jed, area)
879 
880 #ifdef CONSV_HGHT
881 ! Convert tracer moist mixing ratio to mass
882  do iq=1,nwat
883  do k=1,npz
884  do j=js,je
885  do i=is,ie
886  q(i,j,k,iq) = q(i,j,k,iq) * delp(i,j,k)
887  enddo
888  enddo
889  enddo
890  enddo
891 
892  do k=1,npz
893  do j=js,je
894  do i=is,ie
895  ua(i,j,k) = ua(i,j,k) * delp(i,j,k)
896  va(i,j,k) = va(i,j,k) * delp(i,j,k)
897  enddo
898  enddo
899  enddo
900 
901  do j=js,je
902  do i=is,ie
903  pe2(i,1) = ak(1)
904  peln(i,1) = log(pe2(i,1))
905  enddo
906  do k=2,npz+1
907  do i=is,ie
908  pe2(i,k) = pe2(i,k-1) + delp(i,j,k-1)
909  peln(i,k) = log(pe2(i,k))
910  enddo
911  enddo
912  do k=1,npz
913  do i=is,ie
914  pt(i,j,k) = pt(i,j,k)*(peln(i,k+1)-peln(i,k))
915  enddo
916  enddo
917  enddo
918 #endif
919 
920 ! Update ps:
921  do j=js,je
922  do i=is,ie
923  ps(i,j) = ak(1)
924  enddo
925  enddo
926 
927  rdt = dt / (tau_ps/factor + dt)
928  do k=1,npz
929  dbk = rdt*(bk(k+1) - bk(k))
930  do j=js,je
931  do i=is,ie
932  delp(i,j,k) = delp(i,j,k) + dbk*ps_dt(i,j)
933  ps(i,j) = delp(i,j,k) + ps(i,j)
934  enddo
935  enddo
936  enddo
937 
938 #ifdef CONSV_HGHT
939  do j=js,je
940  do k=2,npz+1
941  do i=is,ie
942  pe2(i,k) = pe2(i,k-1) + delp(i,j,k-1)
943  peln(i,k) = log(pe2(i,k))
944  enddo
945  enddo
946  do k=1,npz
947  do i=is,ie
948  pt(i,j,k) = pt(i,j,k)/(peln(i,k+1)-peln(i,k))
949  enddo
950  enddo
951  enddo
952 
953 ! Convert tracer density back to moist mixing ratios
954  do iq=1,nwat
955  do k=1,npz
956  do j=js,je
957  do i=is,ie
958  q(i,j,k,iq) = q(i,j,k,iq) / delp(i,j,k)
959  enddo
960  enddo
961  enddo
962  enddo
963 
964  do k=1,npz
965  do j=js,je
966  do i=is,ie
967  ua(i,j,k) = ua(i,j,k) / delp(i,j,k)
968  va(i,j,k) = va(i,j,k) / delp(i,j,k)
969  enddo
970  enddo
971  enddo
972 #endif
973 
974  end subroutine ps_nudging
975 
976  subroutine ps_bias_correction ( ps_dt, is, ie, js, je, isd, ied, jsd, jed, area )
977  integer, intent(IN) :: is, ie, js, je
978  integer, intent(IN) :: isd, ied, jsd,jed
979  real, intent(inout):: ps_dt(is:ie,js:je)
980  real, intent(IN), dimension(isd:ied,jsd:jed) :: area
981 !
982  real:: esl, total_area
983  real:: bias, psum
984  integer:: i, j
985 
986  total_area = 4.*pi*radius**2
987  esl = 0.01 ! Pascal
988 
989  bias = g0_sum(ps_dt, is, ie, js, je, 1, .true., isd, ied, jsd, jed, area)
990 
991  if ( abs(bias) < esl ) then
992  if(master .and. nudge_debug) write(*,*) 'Very small PS bias=', -bias, ' No bais adjustment'
993  return
994  else
995  if(master .and. nudge_debug) write(*,*) 'Significant PS bias=', -bias
996  endif
997 
998  if ( bias > 0. ) then
999  psum = 0.
1000  do j=js,je
1001  do i=is,ie
1002  if ( ps_dt(i,j) > 0. ) then
1003  psum = psum + area(i,j)
1004  endif
1005  enddo
1006  enddo
1007 
1008  call mp_reduce_sum( psum )
1009  bias = bias * total_area / psum
1010 
1011  do j=js,je
1012  do i=is,ie
1013  if ( ps_dt(i,j) > 0.0 ) then
1014  ps_dt(i,j) = max(0.0, ps_dt(i,j) - bias)
1015  endif
1016  enddo
1017  enddo
1018  else
1019  psum = 0.
1020  do j=js,je
1021  do i=is,ie
1022  if ( ps_dt(i,j) < 0. ) then
1023  psum = psum + area(i,j)
1024  endif
1025  enddo
1026  enddo
1027 
1028  call mp_reduce_sum( psum )
1029  bias = bias * total_area / psum
1030 
1031  do j=js,je
1032  do i=is,ie
1033  if ( ps_dt(i,j) < 0.0 ) then
1034  ps_dt(i,j) = min(0.0, ps_dt(i,j) - bias)
1035  endif
1036  enddo
1037  enddo
1038  endif
1039 
1040  end subroutine ps_bias_correction
1041 
1044  real function g0_sum(p, ifirst, ilast, jfirst, jlast, mode, reproduce, isd, ied, jsd, jed, area)
1045 ! Fast version of global sum; reproduced if result rounded to 4-byte
1046  integer, intent(IN) :: ifirst, ilast
1047  integer, intent(IN) :: jfirst, jlast
1048  integer, intent(IN) :: mode
1049  logical, intent(IN) :: reproduce
1050  real, intent(IN) :: p(ifirst:ilast,jfirst:jlast)
1051  integer, intent(IN) :: isd, ied, jsd, jed
1052  real, intent(IN) :: area(isd:ied,jsd:jed)
1053 
1054  integer :: i,j
1055  real gsum
1056 
1057  gsum = 0.
1058  do j=jfirst,jlast
1059  do i=ifirst,ilast
1060  gsum = gsum + p(i,j)*area(i,j)
1061  enddo
1062  enddo
1063 ! call mpp_sum(gsum) ! does this work?
1064  call mp_reduce_sum(gsum)
1065 
1066  if ( mode==1 ) then
1067  g0_sum = gsum / (4.*pi*radius**2)
1068  else
1069  g0_sum = gsum
1070  endif
1071 
1072  if ( reproduce ) g0_sum = real(g0_sum, 4) ! convert to 4-byte real
1073 
1074  end function g0_sum
1075 
1076 
1077  subroutine compute_slp(isc, iec, jsc, jec, tm, ps, gz, slp)
1078  integer, intent(in):: isc, iec, jsc, jec
1079  real, intent(in), dimension(isc:iec,jsc:jec):: tm, ps, gz
1080 ! tm: virtual temperature required as input
1081  real, intent(out):: slp(isc:iec,jsc:jec)
1082  integer:: i,j
1083 
1084  do j=jsc,jec
1085  do i=isc,iec
1086  slp(i,j) = ps(i,j) * exp( gz(i,j)/(rdgas*(tm(i,j) + 3.25e-3*gz(i,j)/grav)) )
1087  enddo
1088  enddo
1089 
1090  end subroutine compute_slp
1091 
1092 
1093  subroutine get_obs(Time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, pt, nwat, q, u_obs, v_obs, t_obs, q_obs, &
1094  phis, ua, va, u_dt, v_dt, npx, npy, npz, factor, mask, ptop, bd, gridstruct, domain)
1095  type(time_type), intent(in):: Time
1096  integer, intent(in):: npz, nwat, npx, npy
1097  real, intent(in):: zvir, ptop
1098  real, intent(in):: dt, factor
1099  real, intent(in), dimension(npz+1):: ak, bk
1100  type(fv_grid_bounds_type), intent(IN) :: bd
1101  real, intent(in), dimension(isd:ied,jsd:jed):: phis
1102  real, intent(in), dimension(is:ie,js:je):: mask
1103  real, intent(inout), dimension(isd:ied,jsd:jed):: ps
1104  real, intent(inout), dimension(isd:ied,jsd:jed,npz):: delp, pt, ua, va, u_dt, v_dt
1105  real, intent(inout)::q(isd:ied,jsd:jed,npz,*)
1106  real, intent(out), dimension(is:ie,js:je):: ts, ps_obs
1107  real, intent(out), dimension(is:ie,js:je,npz):: u_obs, v_obs, t_obs, q_obs
1108  type(fv_grid_type), intent(IN) :: gridstruct
1109  type(domain2d), intent(INOUT) :: domain
1110 ! local:
1111  real:: tm(is:ie,js:je)
1112  real(KIND=4), allocatable,dimension(:,:,:):: ut, vt, wt
1113  real, allocatable,dimension(:,:,:):: uu, vv
1114  integer :: seconds, days
1115  integer :: i,j,k
1116  real :: alpha, beta
1117 
1118  call get_time (time, seconds, days)
1119 
1120  if ( do_adiabatic_init ) goto 333
1121  seconds = seconds - nint(dt)
1122 
1123 ! Data must be "time_interval" (hr) apart; keep two time levels in memory
1124 
1125  no_obs = .false.
1126  analysis_time = .false.
1127 
1128  if ( mod(seconds, time_interval) == 0 ) then
1129 
1130  if ( nfile == nfile_total ) then
1131  no_obs = .true.
1132 #ifndef DYCORE_SOLO
1133  forecast_mode = .true.
1134 #endif
1135  if(print_end_nudge) then
1136  print_end_nudge = .false.
1137  if (master) write(*,*) '*** L-S nudging Ended at', days, seconds
1138  endif
1139  return ! free-running mode
1140  endif
1141 
1142  ps_dat(:,:,1) = ps_dat(:,:,2)
1143  if ( nudge_winds ) then
1144  u_dat(:,:,:,1) = u_dat(:,:,:,2)
1145  v_dat(:,:,:,1) = v_dat(:,:,:,2)
1146  endif
1147  t_dat(:,:,:,1) = t_dat(:,:,:,2)
1148  q_dat(:,:,:,1) = q_dat(:,:,:,2)
1149 
1150 !---------------
1151 ! Read next data
1152 !---------------
1153  nfile = nfile + 1
1154  call get_ncep_analysis ( ps_dat(:,:,2), u_dat(:,:,:,2), v_dat(:,:,:,2), &
1155  t_dat(:,:,:,2), q_dat(:,:,:,2), zvir, &
1156  ts, nfile, file_names(nfile), bd )
1157  time_nudge = dt
1158  else
1159  time_nudge = time_nudge + dt
1160  endif
1161 
1162 !--------------------
1163 ! Time interpolation:
1164 !--------------------
1165 
1166  beta = time_nudge / real(time_interval)
1167 
1168  if ( beta < 0. .or. beta > (1.+1.e-7) ) then
1169  if (master) write(*,*) 'Nudging: beta=', beta
1170  call mpp_error(fatal,'==> Error from get_obs:data out of range')
1171  endif
1172 
1173 333 continue
1174 
1175  if ( do_adiabatic_init ) then
1176  beta = 1.; alpha = 0.
1177  if( nudge_debug .and. master) write(*,*) 'Doing final adiabatic initialization/nudging'
1178  else
1179  alpha = 1. - beta
1180  if( abs(alpha)<1.e-7 ) analysis_time = .true.
1181  endif
1182 
1183 ! Warning: ps_data are not adjusted for the differences in terrain yet
1184  ps_obs(:,:) = alpha*ps_dat(:,:,1) + beta*ps_dat(:,:,2)
1185 
1186 !---------------------------------
1187 !*** nudge & update ps & delp here
1188 !---------------------------------
1189  if (nudge_ps) then
1190 
1191  allocate ( wt(is:ie,js:je,km) )
1192  wt(:,:,:) = alpha*t_dat(:,:,:,1) + beta*t_dat(:,:,:,2)
1193 ! Needs gz3 for ps_nudging
1194  call get_int_hght(npz, ak, bk, ps(is:ie,js:je), delp, ps_obs(is:ie,js:je), wt)
1195  do j=js,je
1196  do i=is,ie
1197  tm(i,j) = wt(i,j,km)
1198  enddo
1199  enddo
1200  deallocate ( wt )
1201 
1202  allocate ( uu(isd:ied,jsd:jed,npz) )
1203  allocate ( vv(isd:ied,jsd:jed,npz) )
1204  uu = ua
1205  vv = va
1206  call ps_nudging(dt, factor, npz, ak, bk, ps_obs, mask, tm, ps, phis, delp, uu, vv, pt, nwat, q, bd, npx, npy, gridstruct, domain)
1207  do k=1,npz
1208  do j=js,je
1209  do i=is,ie
1210  u_dt(i,j,k) = u_dt(i,j,k) + (uu(i,j,k) - ua(i,j,k)) / dt
1211  v_dt(i,j,k) = v_dt(i,j,k) + (vv(i,j,k) - va(i,j,k)) / dt
1212  enddo
1213  enddo
1214  enddo
1215  deallocate (uu )
1216  deallocate (vv )
1217  endif
1218 
1219  allocate ( ut(is:ie,js:je,npz) )
1220  allocate ( vt(is:ie,js:je,npz) )
1221 
1222  if ( nudge_winds ) then
1223 
1224  call remap_uv(npz, ak, bk, ps(is:ie,js:je), delp, ut, vt, &
1225  km, ps_dat(is:ie,js:je,1), u_dat(:,:,:,1), v_dat(:,:,:,1), ptop )
1226 
1227  u_obs(:,:,:) = alpha*ut(:,:,:)
1228  v_obs(:,:,:) = alpha*vt(:,:,:)
1229 
1230  call remap_uv(npz, ak, bk, ps(is:ie,js:je), delp, ut, vt, &
1231  km, ps_dat(is:ie,js:je,2), u_dat(:,:,:,2), v_dat(:,:,:,2), ptop )
1232 
1233  u_obs(:,:,:) = u_obs(:,:,:) + beta*ut(:,:,:)
1234  v_obs(:,:,:) = v_obs(:,:,:) + beta*vt(:,:,:)
1235  endif
1236 
1237  call remap_tq(npz, ak, bk, ps(is:ie,js:je), delp, ut, vt, &
1238  km, ps_dat(is:ie,js:je,1), t_dat(:,:,:,1), q_dat(:,:,:,1), zvir, ptop)
1239 
1240  t_obs(:,:,:) = alpha*ut(:,:,:)
1241  q_obs(:,:,:) = alpha*vt(:,:,:)
1242 
1243  call remap_tq(npz, ak, bk, ps(is:ie,js:je), delp, ut, vt, &
1244  km, ps_dat(is:ie,js:je,2), t_dat(:,:,:,2), q_dat(:,:,:,2), zvir, ptop)
1245 
1246  t_obs(:,:,:) = t_obs(:,:,:) + beta*ut(:,:,:)
1247  q_obs(:,:,:) = q_obs(:,:,:) + beta*vt(:,:,:)
1248 
1249  deallocate ( ut )
1250  deallocate ( vt )
1251 
1252  end subroutine get_obs
1253 
1256  subroutine fv_nwp_nudge_init(time, axes, npz, zvir, ak, bk, ts, phis, gridstruct, ks, npx, neststruct, bd)
1257  character(len=17) :: mod_name = 'fv_nudge'
1258  type(time_type), intent(in):: time
1259  integer, intent(in):: axes(4)
1260  integer, intent(in):: npz
1261  real, intent(in):: zvir
1262  type(fv_grid_bounds_type), intent(IN) :: bd
1263  real, intent(in), dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: phis
1264  real, intent(in), dimension(npz+1):: ak, bk
1265  real, intent(out), dimension(bd%is:bd%ie,bd%js:bd%je):: ts
1266  type(fv_grid_type), target :: gridstruct
1267  integer, intent(in) :: ks, npx
1268  type(fv_nest_type) :: neststruct
1269 
1270  real :: missing_value = -1.e10
1271  logical found
1272  integer tsize(4)
1273  integer :: i, j, j1, f_unit, unit, io, ierr, nt, k
1274  integer :: ncid
1275 
1276 ! ---> h1g, read NCEP analysis data list, 2012-10-22
1277  integer :: input_fname_unit
1278  character(len=128) :: fname_tmp
1279 ! <--- h1g, 2012-10-22
1280 
1281  real, pointer, dimension(:,:,:) :: agrid
1282 
1283  is = bd%is
1284  ie = bd%ie
1285  js = bd%js
1286  je = bd%je
1287 
1288  isd = bd%isd
1289  ied = bd%ied
1290  jsd = bd%jsd
1291  jed = bd%jed
1292 
1293 
1294  agrid => gridstruct%agrid
1295 
1296  master = is_master()
1297  do_adiabatic_init = .false.
1298  deg2rad = pi/180.
1299  rad2deg = 180./pi
1300 
1301  if (neststruct%nested) then
1302 !!! Assumes no grid stretching
1303  grid_size = 1.e7/(neststruct%npx_global*neststruct%refinement_of_global)
1304  else
1305  grid_size = 1.e7/real(npx-1) ! mean grid size
1306  endif
1307 
1308  do nt=1,nfile_max
1309  file_names(nt) = "No_File_specified"
1310  enddo
1311 
1312  track_file_name = "No_File_specified"
1313 
1314  if( file_exist( 'input.nml' ) ) then
1315  unit = open_namelist_file()
1316  io = 1
1317  do while ( io .ne. 0 )
1318  read( unit, nml = fv_nwp_nudge_nml, iostat = io, end = 10 )
1319  ierr = check_nml_error(io,'fv_nwp_nudge_nml')
1320  end do
1321 10 call close_file ( unit )
1322  end if
1323  call write_version_number ( 'FV_NUDGE_MOD', version )
1324  if ( master ) then
1325  f_unit=stdlog()
1326  write( f_unit, nml = fv_nwp_nudge_nml )
1327  write(*,*) 'NWP nudging initialized.'
1328  endif
1329 ! ---> h1g, specify the NCEP analysis data for nudging, 2012-10-22
1330  if ( trim(input_fname_list) .ne. "" ) then
1331  input_fname_unit = get_unit()
1332  open( input_fname_unit, file = input_fname_list )
1333  nt = 0
1334  io = 0
1335 
1336  do while ( io .eq. 0 )
1337  read ( input_fname_unit, '(a)', iostat = io, end = 101 ) fname_tmp
1338  if( trim(fname_tmp) .ne. "" ) then ! escape any empty record
1339  if ( trim(fname_tmp) == trim(analysis_file_last) ) then
1340  nt = nt + 1
1341  file_names(nt) = 'INPUT/'//trim(fname_tmp)
1342  if(master .and. nudge_debug) write(*,*) 'From NCEP file list, last file: ', nt, file_names(nt)
1343  nt = 0
1344  goto 101 ! read last analysis data and then close file
1345  endif ! trim(fname_tmp) == trim(analysis_file_last)
1346 
1347  if ( trim(analysis_file_first) == "" ) then
1348  nt = nt + 1
1349  file_names(nt) = 'INPUT/'//trim(fname_tmp)
1350  if(master .and. nudge_debug) then
1351  if( nt .eq. 1 ) then
1352  write(*,*) 'From NCEP file list, first file: ', nt, file_names(nt),trim(analysis_file_first)
1353  else
1354  write(*,*) 'From NCEP file list: ', nt, file_names(nt)
1355  endif ! nt .eq. 1
1356  endif ! master .and. nudge_debug
1357  else
1358  if ( trim(fname_tmp) == trim(analysis_file_first) .or. nt > 0 ) then
1359  nt = nt + 1
1360  file_names(nt) = 'INPUT/'//trim(fname_tmp)
1361  if(master .and. nudge_debug) then
1362  if( nt .eq. 1 ) then
1363  write(*,*) 'From NCEP file list, first file: ', nt, file_names(nt),trim(analysis_file_first)
1364  else
1365  write(*,*) 'From NCEP file list: ', nt, file_names(nt)
1366  endif ! nt .eq. 1
1367  endif ! master .and. nudge_debug
1368  endif ! trim(fname_tmp) == trim(analysis_file_first) .or. nt > 0
1369  endif ! trim(analysis_file_first) == ""
1370  endif ! trim(fname_tmp) .ne. ""
1371  end do ! io .eq. 0
1372 101 close( input_fname_unit )
1373  endif
1374 ! <--- h1g, 2012-10-22
1375 
1376  do nt=1,nfile_max
1377  if ( file_names(nt) == "No_File_specified" ) then
1378  nfile_total = nt - 1
1379  if(master) write(*,*) 'Total of NCEP files specified=', nfile_total
1380  exit
1381  endif
1382  enddo
1383 
1384  id_ht_err = register_diag_field( mod_name, 'ht_error', axes(1:2), time, &
1385  'height_error', 'DEG K', missing_value=missing_value )
1386 
1387 
1388 ! Initialize remapping coefficients:
1389 
1390 ! call field_size(file_names(1), 'T', tsize, field_found=found)
1391 ! if ( found ) then
1392 ! im = tsize(1); jm = tsize(2); km = tsize(3)
1393 ! if(master) write(*,*) 'NCEP analysis dimensions:', tsize
1394 ! else
1395 ! call mpp_error(FATAL,'==> Error from get_ncep_analysis: T field not found')
1396 ! endif
1397  call open_ncfile( file_names(1), ncid ) ! open the file
1398  call get_ncdim1( ncid, 'lon', im )
1399  call get_ncdim1( ncid, 'lat', jm )
1400  call get_ncdim1( ncid, 'lev', km )
1401  if(master) write(*,*) 'NCEP analysis dimensions:', im, jm, km
1402 
1403  allocate ( s2c(is:ie,js:je,4) )
1404  allocate ( id1(is:ie,js:je) )
1405  allocate ( id2(is:ie,js:je) )
1406  allocate ( jdc(is:ie,js:je) )
1407 
1408  allocate ( lon(im) )
1409  allocate ( lat(jm) )
1410 
1411  call _get_var1 (ncid, 'lon', im, lon )
1412  call _get_var1 (ncid, 'lat', jm, lat )
1413 
1414 ! Convert to radian
1415  do i=1,im
1416  lon(i) = lon(i) * deg2rad ! lon(1) = 0.
1417  enddo
1418  do j=1,jm
1419  lat(j) = lat(j) * deg2rad
1420  enddo
1421 
1422  allocate ( ak0(km+1) )
1423  allocate ( bk0(km+1) )
1424 
1425  call _get_var1 (ncid, 'hyai', km+1, ak0, found )
1426  if ( .not. found ) ak0(:) = 0.
1427 
1428  call _get_var1 (ncid, 'hybi', km+1, bk0 )
1429  call close_ncfile( ncid )
1430 
1431 ! Note: definition of NCEP hybrid is p(k) = a(k)*1.E5 + b(k)*ps
1432  ak0(:) = ak0(:) * 1.e5
1433 ! Limiter to prevent NAN at top during remapping
1434  if ( bk0(1) < 1.e-9 ) ak0(1) = max(1.e-9, ak0(1))
1435 
1436  if ( master ) then
1437  do k=1,npz
1438  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)
1439  enddo
1440  endif
1441 
1442  if ( k_breed==0 ) k_breed = max(1, ks)
1443 
1444  call slp_obs_init
1445 
1446 !-----------------------------------------------------------
1447 ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
1448 !-----------------------------------------------------------
1449  call remap_coef(agrid)
1450 
1451 ! Find bounding latitudes:
1452  jbeg = jm-1; jend = 2
1453  do j=js,je
1454  do i=is,ie
1455  j1 = jdc(i,j)
1456  jbeg = min(jbeg, j1)
1457  jend = max(jend, j1+1)
1458  enddo
1459  enddo
1460 
1461  allocate ( gz0(is:ie,js:je) )
1462  allocate ( gz3(is:ie,js:je,km+1) )
1463  allocate (ps_dat(is:ie,js:je,2) )
1464  allocate ( u_dat(is:ie,js:je,km,2) )
1465  allocate ( v_dat(is:ie,js:je,km,2) )
1466  allocate ( t_dat(is:ie,js:je,km,2) )
1467  allocate ( q_dat(is:ie,js:je,km,2) )
1468 
1469 
1470 ! Get first dataset
1471  nt = 2
1472  nfile = 1
1473  call get_ncep_analysis ( ps_dat(:,:,nt), u_dat(:,:,:,nt), v_dat(:,:,:,nt), &
1474  t_dat(:,:,:,nt), q_dat(:,:,:,nt), zvir, &
1475  ts, nfile, file_names(nfile), bd )
1476 
1477 
1478  module_is_initialized = .true.
1479 
1480  nullify(agrid)
1481 
1482  end subroutine fv_nwp_nudge_init
1483 
1484 
1485  subroutine get_ncep_analysis ( ps, u, v, t, q, zvir, ts, nfile, fname, bd )
1486  real, intent(in):: zvir
1487  character(len=128), intent(in):: fname
1488  integer, intent(inout):: nfile
1489 !
1490  type(fv_grid_bounds_type), intent(IN) :: bd
1491  real, intent(out), dimension(is:ie,js:je):: ts, ps
1492  real(KIND=4), intent(out), dimension(is:ie,js:je,km):: u, v, t, q
1493 ! local:
1494  real(kind=4), allocatable:: wk0(:,:), wk1(:,:), wk2(:,:), wk3(:,:,:)
1495  real tmean
1496  integer:: i, j, k, npt
1497  integer:: i1, i2, j1, ncid
1498  logical found
1499  logical:: read_ts = .true.
1500  logical:: land_ts = .false.
1501 
1502  integer:: status, var3id ! h1g, 2016-08-10
1503 #include <netcdf.inc>
1504 
1505 
1506  if( .not. file_exist(fname) ) then
1507  call mpp_error(fatal,'==> Error from get_ncep_analysis: file not found')
1508  else
1509  call open_ncfile( fname, ncid ) ! open the file
1510  if(master) write(*,*) 'Reading NCEP anlysis file:', fname
1511  endif
1512 
1513  if ( read_ts ) then ! read skin temperature; could be used for SST
1514  allocate ( wk1(im,jm) )
1515 
1516  call get_var3_r4( ncid, 'TS', 1,im, 1,jm, 1,1, wk1 )
1517 ! if ( master ) write(*,*) 'Done reading NCEP TS data'
1518 
1519  if ( .not. land_ts ) then
1520  allocate ( wk0(im,jm) )
1521 ! Read NCEP ORO (1; land; 0: ocean; 2: sea_ice)
1522 
1523 ! ---> h1g, read either 'ORO' or 'LAND', 2016-08-10
1524  status = nf_inq_varid(ncid, 'ORO', var3id)
1525  if (status .eq. nf_noerr) then
1526  call get_var3_r4( ncid, 'ORO', 1,im, 1,jm, 1,1, wk0 )
1527 
1528  else !there is no 'ORO'
1529  status = nf_inq_varid(ncid, 'LAND', var3id)
1530  if (status .eq. nf_noerr) then
1531  call get_var3_r4( ncid, 'LAND', 1,im, 1,jm, 1,1, wk0 )
1532  else
1533  call mpp_error(fatal,'Neither ORO nor LAND exists in re-analysis data')
1534  endif
1535  endif
1536 ! <--- h1g, 2016-08-10
1537 
1538  do j=1,jm
1539  tmean = 0.
1540  npt = 0
1541  do i=1,im
1542  if( abs(wk0(i,j)-1.) > 0.99 ) then ! non-land points
1543  tmean = tmean + wk1(i,j)
1544  npt = npt + 1
1545  endif
1546  enddo
1547 !-------------------------------------------------------
1548 ! Replace TS over interior land with zonal mean SST/Ice
1549 !-------------------------------------------------------
1550  if ( npt /= 0 ) then
1551  tmean= tmean / real(npt)
1552  do i=1,im
1553  if( abs(wk0(i,j)-1.) <= 0.99 ) then ! land points
1554  if ( i==1 ) then
1555  i1 = im; i2 = 2
1556  elseif ( i==im ) then
1557  i1 = im-1; i2 = 1
1558  else
1559  i1 = i-1; i2 = i+1
1560  endif
1561  if ( abs(wk0(i2,j)-1.)>0.99 ) then ! east side has priority
1562  wk1(i,j) = wk1(i2,j)
1563  elseif ( abs(wk0(i1,j)-1.)>0.99 ) then ! west side
1564  wk1(i,j) = wk1(i1,j)
1565  else
1566  wk1(i,j) = tmean
1567  endif
1568  endif
1569  enddo
1570  endif
1571  enddo
1572  deallocate ( wk0 )
1573  endif ! land_ts
1574 
1575  do j=js,je
1576  do i=is,ie
1577  i1 = id1(i,j)
1578  i2 = id2(i,j)
1579  j1 = jdc(i,j)
1580  ts(i,j) = s2c(i,j,1)*wk1(i1,j1 ) + s2c(i,j,2)*wk1(i2,j1 ) + &
1581  s2c(i,j,3)*wk1(i2,j1+1) + s2c(i,j,4)*wk1(i1,j1+1)
1582  enddo
1583  enddo
1584  call prt_maxmin('SST_model', ts, is, ie, js, je, 0, 1, 1.)
1585 
1586 #ifndef DYCORE_SOLO
1587 ! Perform interp to FMS SST format/grid
1588  call ncep2fms( wk1 )
1589  if(master) call pmaxmin( 'SST_ncep', sst_ncep, i_sst, j_sst, 1.)
1590 ! if(nfile/=1 .and. master) call pmaxmin( 'SST_anom', sst_anom, i_sst, j_sst, 1.)
1591 #endif
1592  deallocate ( wk1 )
1593  if (master) write(*,*) 'Done processing NCEP SST'
1594 
1595  endif ! read_ts
1596 
1597 !----------------------------------
1598 ! remap surface pressure and height:
1599 !----------------------------------
1600 
1601  allocate ( wk2(im,jbeg:jend) )
1602  call get_var3_r4( ncid, 'PS', 1,im, jbeg,jend, 1,1, wk2 )
1603 
1604  do j=js,je
1605  do i=is,ie
1606  i1 = id1(i,j)
1607  i2 = id2(i,j)
1608  j1 = jdc(i,j)
1609  ps(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
1610  s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
1611  enddo
1612  enddo
1613 
1614 
1615 ! ---> h1g, read either 'PHIS' or 'PHI', 2016-08-10
1616  status = nf_inq_varid(ncid, 'PHIS', var3id)
1617  if (status .eq. nf_noerr) then
1618  call get_var3_r4( ncid, 'PHIS', 1,im, jbeg,jend, 1,1, wk2 )
1619 
1620  else !there is no 'PHIS'
1621  status = nf_inq_varid(ncid, 'PHI', var3id)
1622  if (status .eq. nf_noerr) then
1623  call get_var3_r4( ncid, 'PHI', 1,im, jbeg,jend, 1,1, wk2 )
1624  wk2 = wk2 * grav ! convert unit from geopotential meter (m) to geopotential height (m2/s2)
1625  else
1626  call mpp_error(fatal,'Neither PHIS nor PHI exists in re-analysis data')
1627  endif
1628  endif
1629 ! <--- h1g, 2016-08-10
1630 
1631 
1632  do j=js,je
1633  do i=is,ie
1634  i1 = id1(i,j)
1635  i2 = id2(i,j)
1636  j1 = jdc(i,j)
1637  gz0(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
1638  s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
1639  enddo
1640  enddo
1641  call prt_maxmin('ZS_ncep', gz0, is, ie, js, je, 0, 1, 1./grav)
1642  deallocate ( wk2 )
1643 
1644 
1645  allocate ( wk3(1:im, jbeg:jend, 1:km) )
1646 
1647 ! Winds:
1648  if ( nudge_winds ) then
1649 
1650  call get_var3_r4( ncid, 'U', 1,im, jbeg,jend, 1,km, wk3 )
1651 
1652  do k=1,km
1653  do j=js,je
1654  do i=is,ie
1655  i1 = id1(i,j)
1656  i2 = id2(i,j)
1657  j1 = jdc(i,j)
1658  u(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1659  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1660  enddo
1661  enddo
1662  enddo
1663 
1664  call get_var3_r4( ncid, 'V', 1,im, jbeg,jend, 1,km, wk3 )
1665 
1666  do k=1,km
1667  do j=js,je
1668  do i=is,ie
1669  i1 = id1(i,j)
1670  i2 = id2(i,j)
1671  j1 = jdc(i,j)
1672  v(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1673  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1674  enddo
1675  enddo
1676  enddo
1677 
1678  endif
1679 
1680 ! Read in tracers: only sphum at this point
1681  call get_var3_r4( ncid, 'Q', 1,im, jbeg,jend, 1,km , wk3 )
1682 
1683  do k=1,km
1684  do j=js,je
1685  do i=is,ie
1686  i1 = id1(i,j)
1687  i2 = id2(i,j)
1688  j1 = jdc(i,j)
1689  q(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1690  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1691  enddo
1692  enddo
1693  enddo
1694 
1695  call get_var3_r4( ncid, 'T', 1,im, jbeg,jend, 1,km , wk3 )
1696  call close_ncfile ( ncid )
1697 
1698  do k=1,km
1699  do j=js,je
1700  do i=is,ie
1701  i1 = id1(i,j)
1702  i2 = id2(i,j)
1703  j1 = jdc(i,j)
1704  t(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1705  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1706  enddo
1707  enddo
1708  enddo
1709 
1710  if ( .not. t_is_tv ) then
1711  do k=1,km
1712  do j=js,je
1713  do i=is,ie
1714 ! The field T in Larry H.'s post processing of NCEP analysis is actually virtual temperature
1715 ! before Dec 1, 2005
1716 ! Convert t to virtual temperature:
1717  t(i,j,k) = t(i,j,k)*(1.+zvir*q(i,j,k))
1718  enddo
1719  enddo
1720  enddo
1721  endif
1722 
1723 ! endif
1724 
1725  deallocate ( wk3 )
1726 
1727 ! nfile = nfile + 1
1728 
1729  end subroutine get_ncep_analysis
1730 
1731 
1732  subroutine remap_coef(agrid)
1734  real, intent(IN), dimension(isd:ied,jsd:jed,2) :: agrid
1735 
1736 ! local:
1737  real :: rdlon(im)
1738  real :: rdlat(jm)
1739  real:: a1, b1
1740  integer i,j, i1, i2, jc, i0, j0
1741 
1742  do i=1,im-1
1743  rdlon(i) = 1. / (lon(i+1) - lon(i))
1744  enddo
1745  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
1746 
1747  do j=1,jm-1
1748  rdlat(j) = 1. / (lat(j+1) - lat(j))
1749  enddo
1750 
1751 ! * Interpolate to cubed sphere cell center
1752  do 5000 j=js,je
1753 
1754  do i=is,ie
1755 
1756  if ( agrid(i,j,1)>lon(im) ) then
1757  i1 = im; i2 = 1
1758  a1 = (agrid(i,j,1)-lon(im)) * rdlon(im)
1759  elseif ( agrid(i,j,1)<lon(1) ) then
1760  i1 = im; i2 = 1
1761  a1 = (agrid(i,j,1)+2.*pi-lon(im)) * rdlon(im)
1762  else
1763  do i0=1,im-1
1764  if ( agrid(i,j,1)>=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then
1765  i1 = i0; i2 = i0+1
1766  a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0)
1767  go to 111
1768  endif
1769  enddo
1770  endif
1771 111 continue
1772 
1773  if ( agrid(i,j,2)<lat(1) ) then
1774  jc = 1
1775  b1 = 0.
1776  elseif ( agrid(i,j,2)>lat(jm) ) then
1777  jc = jm-1
1778  b1 = 1.
1779  else
1780  do j0=1,jm-1
1781  if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then
1782  jc = j0
1783  b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc)
1784  go to 222
1785  endif
1786  enddo
1787  endif
1788 222 continue
1789 
1790  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
1791  write(*,*) 'gid=', mpp_pe(), i,j,a1, b1
1792  endif
1793 
1794  s2c(i,j,1) = (1.-a1) * (1.-b1)
1795  s2c(i,j,2) = a1 * (1.-b1)
1796  s2c(i,j,3) = a1 * b1
1797  s2c(i,j,4) = (1.-a1) * b1
1798  id1(i,j) = i1
1799  id2(i,j) = i2
1800  jdc(i,j) = jc
1801  enddo !i-loop
1802 5000 continue ! j-loop
1803 
1804  end subroutine remap_coef
1805 
1806 
1807 #ifndef DYCORE_SOLO
1808  subroutine ncep2fms( sst )
1809  real(kind=4), intent(in):: sst(im,jm)
1810 ! local:
1811  real :: rdlon(im)
1812  real :: rdlat(jm)
1813  real:: a1, b1
1814  real:: delx, dely
1815  real:: xc, yc ! "data" location
1816  real:: c1, c2, c3, c4
1817  integer i,j, i1, i2, jc, i0, j0, it, jt
1818 
1819  do i=1,im-1
1820  rdlon(i) = 1. / (lon(i+1) - lon(i))
1821  enddo
1822  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
1823 
1824  do j=1,jm-1
1825  rdlat(j) = 1. / (lat(j+1) - lat(j))
1826  enddo
1827 
1828 ! * Interpolate to "FMS" 1x1 SST data grid
1829 ! lon: 0.5, 1.5, ..., 359.5
1830 ! lat: -89.5, -88.5, ... , 88.5, 89.5
1831 
1832  delx = 360./real(i_sst)
1833  dely = 180./real(j_sst)
1834 
1835  jt = 1
1836  do 5000 j=1,j_sst
1837 
1838  yc = (-90. + dely * (0.5+real(j-1))) * deg2rad
1839  if ( yc<lat(1) ) then
1840  jc = 1
1841  b1 = 0.
1842  elseif ( yc>lat(jm) ) then
1843  jc = jm-1
1844  b1 = 1.
1845  else
1846  do j0=jt,jm-1
1847  if ( yc>=lat(j0) .and. yc<=lat(j0+1) ) then
1848  jc = j0
1849  jt = j0
1850  b1 = (yc-lat(jc)) * rdlat(jc)
1851  go to 222
1852  endif
1853  enddo
1854  endif
1855 222 continue
1856  it = 1
1857 
1858  do i=1,i_sst
1859  xc = delx * (0.5+real(i-1)) * deg2rad
1860  if ( xc>lon(im) ) then
1861  i1 = im; i2 = 1
1862  a1 = (xc-lon(im)) * rdlon(im)
1863  elseif ( xc<lon(1) ) then
1864  i1 = im; i2 = 1
1865  a1 = (xc+2.*pi-lon(im)) * rdlon(im)
1866  else
1867  do i0=it,im-1
1868  if ( xc>=lon(i0) .and. xc<=lon(i0+1) ) then
1869  i1 = i0; i2 = i0+1
1870  it = i0
1871  a1 = (xc-lon(i1)) * rdlon(i0)
1872  go to 111
1873  endif
1874  enddo
1875  endif
1876 111 continue
1877 
1878 ! if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
1879 ! write(*,*) 'gid=', mpp_pe(), i,j,a1, b1
1880 ! endif
1881  c1 = (1.-a1) * (1.-b1)
1882  c2 = a1 * (1.-b1)
1883  c3 = a1 * b1
1884  c4 = (1.-a1) * b1
1885 ! Interpolated surface pressure
1886  sst_ncep(i,j) = c1*sst(i1,jc ) + c2*sst(i2,jc ) + &
1887  c3*sst(i2,jc+1) + c4*sst(i1,jc+1)
1888  enddo !i-loop
1889 5000 continue ! j-loop
1890 
1891  end subroutine ncep2fms
1892 #endif
1893 
1894  subroutine get_int_hght(npz, ak, bk, ps, delp, ps0, tv)
1895  integer, intent(in):: npz
1896  real, intent(in):: ak(npz+1), bk(npz+1)
1897  real, intent(in), dimension(is:ie,js:je):: ps, ps0
1898  real, intent(in), dimension(isd:ied,jsd:jed,npz):: delp
1899  real(KIND=4), intent(in), dimension(is:ie,js:je,km):: tv ! virtual temperature
1900 ! local:
1901  real, dimension(is:ie,km+1):: pn0
1902  integer i,j,k
1903 
1904 !$OMP parallel do default(none) shared(is,ie,js,je,km,ak0,bk0,ps0,gz3,gz0,tv) &
1905 !$OMP private(pn0)
1906  do 5000 j=js,je
1907 
1908  do k=1,km+1
1909  do i=is,ie
1910  pn0(i,k) = log( ak0(k) + bk0(k)*ps0(i,j) )
1911  enddo
1912  enddo
1913  do i=is,ie
1914  gz3(i,j,km+1) = gz0(i,j) ! Data Surface geopotential
1915  enddo
1916 
1917  do k=km,1,-1
1918  do i=is,ie
1919  gz3(i,j,k) = gz3(i,j,k+1) + rdgas*tv(i,j,k)*(pn0(i,k+1)-pn0(i,k))
1920  enddo
1921  enddo
1922 
1923 5000 continue
1924 
1925  end subroutine get_int_hght
1926 
1927 
1928 
1929  subroutine remap_tq( npz, ak, bk, ps, delp, t, q, &
1930  kmd, ps0, ta, qa, zvir, ptop)
1931  integer, intent(in):: npz, kmd
1932  real, intent(in):: zvir, ptop
1933  real, intent(in):: ak(npz+1), bk(npz+1)
1934  real, intent(in), dimension(is:ie,js:je):: ps0
1935  real, intent(inout), dimension(is:ie,js:je):: ps
1936  real, intent(in), dimension(isd:ied,jsd:jed,npz):: delp
1937  real(KIND=4), intent(in), dimension(is:ie,js:je,kmd):: ta, qa
1938  real(KIND=4), intent(out), dimension(is:ie,js:je,npz):: t
1939 ! on input "ta" is virtual temperature
1940 ! on output "t" is virtual temperature
1941  real(KIND=4), intent(out), dimension(is:ie,js:je,npz):: q
1942 ! local:
1943  real, dimension(is:ie,kmd):: tp, qp
1944  real, dimension(is:ie,kmd+1):: pe0, pn0
1945  real, dimension(is:ie,npz):: qn1
1946  real, dimension(is:ie,npz+1):: pe1, pn1
1947  integer i,j,k
1948 
1949  do 5000 j=js,je
1950 
1951  do k=1,kmd+1
1952  do i=is,ie
1953  pe0(i,k) = ak0(k) + bk0(k)*ps0(i,j)
1954  pn0(i,k) = log(pe0(i,k))
1955  enddo
1956  enddo
1957 !------
1958 ! Model
1959 !------
1960  do i=is,ie
1961  pe1(i,1) = ak(1)
1962  enddo
1963  do k=1,npz
1964  do i=is,ie
1965  pe1(i,k+1) = pe1(i,k) + delp(i,j,k)
1966  enddo
1967  enddo
1968  do i=is,ie
1969  ps(i,j) = pe1(i,npz+1)
1970  enddo
1971  do k=1,npz+1
1972  do i=is,ie
1973  pn1(i,k) = log(pe1(i,k))
1974  enddo
1975  enddo
1976 
1977  if ( nudge_q ) then
1978  do k=1,kmd
1979  do i=is,ie
1980  qp(i,k) = qa(i,j,k)
1981  enddo
1982  enddo
1983  call mappm(kmd, pe0, qp, npz, pe1, qn1, is,ie, 0, kord_data, ptop)
1984  do k=1,npz
1985  do i=is,ie
1986  q(i,j,k) = qn1(i,k)
1987  enddo
1988  enddo
1989  endif
1990 
1991  do k=1,kmd
1992  do i=is,ie
1993  tp(i,k) = ta(i,j,k)
1994  enddo
1995  enddo
1996  call mappm(kmd, pn0, tp, npz, pn1, qn1, is,ie, 1, kord_data, ptop)
1997 
1998  do k=1,npz
1999  do i=is,ie
2000  t(i,j,k) = qn1(i,k)
2001  enddo
2002  enddo
2003 
2004 5000 continue
2005 
2006  end subroutine remap_tq
2007 
2008 
2009  subroutine remap_uv(npz, ak, bk, ps, delp, u, v, kmd, ps0, u0, v0, ptop)
2010  integer, intent(in):: npz
2011  real, intent(IN):: ptop
2012  real, intent(in):: ak(npz+1), bk(npz+1)
2013  real, intent(inout):: ps(is:ie,js:je)
2014  real, intent(in), dimension(isd:ied,jsd:jed,npz):: delp
2015  real(KIND=4), intent(inout), dimension(is:ie,js:je,npz):: u, v
2016 !
2017  integer, intent(in):: kmd
2018  real, intent(in):: ps0(is:ie,js:je)
2019  real(KIND=4), intent(in), dimension(is:ie,js:je,kmd):: u0, v0
2020 !
2021 ! local:
2022  real, dimension(is:ie,kmd+1):: pe0
2023  real, dimension(is:ie,npz+1):: pe1
2024  real, dimension(is:ie,kmd):: qt
2025  real, dimension(is:ie,npz):: qn1
2026  integer i,j,k
2027 
2028  do 5000 j=js,je
2029 !------
2030 ! Data
2031 !------
2032  do k=1,kmd+1
2033  do i=is,ie
2034  pe0(i,k) = ak0(k) + bk0(k)*ps0(i,j)
2035  enddo
2036  enddo
2037 !------
2038 ! Model
2039 !------
2040  do i=is,ie
2041  pe1(i,1) = ak(1)
2042  enddo
2043  do k=1,npz
2044  do i=is,ie
2045  pe1(i,k+1) = pe1(i,k) + delp(i,j,k)
2046  enddo
2047  enddo
2048  do i=is,ie
2049  ps(i,j) = pe1(i,npz+1)
2050  enddo
2051 !------
2052 ! map u
2053 !------
2054  do k=1,kmd
2055  do i=is,ie
2056  qt(i,k) = u0(i,j,k)
2057  enddo
2058  enddo
2059  call mappm(kmd, pe0, qt, npz, pe1, qn1, is,ie, -1, kord_data, ptop)
2060  do k=1,npz
2061  do i=is,ie
2062  u(i,j,k) = qn1(i,k)
2063  enddo
2064  enddo
2065 !------
2066 ! map v
2067 !------
2068  do k=1,kmd
2069  do i=is,ie
2070  qt(i,k) = v0(i,j,k)
2071  enddo
2072  enddo
2073  call mappm(kmd, pe0, qt, npz, pe1, qn1, is,ie, -1, kord_data, ptop)
2074  do k=1,npz
2075  do i=is,ie
2076  v(i,j,k) = qn1(i,k)
2077  enddo
2078  enddo
2079 5000 continue
2080 
2081  end subroutine remap_uv
2082 
2084  subroutine fv_nwp_nudge_end
2085  deallocate ( ps_dat )
2086  deallocate ( t_dat )
2087  deallocate ( q_dat )
2088 
2089  if ( nudge_winds ) then
2090  deallocate ( u_dat )
2091  deallocate ( v_dat )
2092  endif
2093 
2094  deallocate ( s2c )
2095  deallocate ( id1 )
2096  deallocate ( id2 )
2097  deallocate ( jdc )
2098 
2099  deallocate ( ak0 )
2100  deallocate ( bk0 )
2101  deallocate ( lat )
2102  deallocate ( lon )
2103 
2104  deallocate ( gz3 )
2105  deallocate ( gz0 )
2106 
2107  end subroutine fv_nwp_nudge_end
2108 
2109 
2110  subroutine get_tc_mask(time, mask, agrid)
2111  real :: slp_mask = 100900.
2112 ! Input
2113  type(time_type), intent(in):: time
2114  real, intent(inout):: mask(is:ie,js:je)
2115  real(kind=R_GRID), intent(IN), dimension(isd:ied,jsd:jed,2) :: agrid
2116 ! local
2117  real(kind=R_GRID):: pos(2)
2118  real:: slp_o
2119  real:: w10_o
2120  real:: r_vor, p_vor
2121  real:: dist
2122  integer n, i, j
2123 
2124  do 5000 n=1,nstorms ! loop through all storms
2125 !----------------------------------------
2126 ! Obtain slp observation
2127 !----------------------------------------
2128  call get_slp_obs(time, nobs_tc(n), x_obs(1,n), y_obs(1,n), wind_obs(1,n), mslp_obs(1,n), mslp_out(1,n), rad_out(1,n), &
2129  time_tc(1,n), pos(1), pos(2), w10_o, slp_o, r_vor, p_vor)
2130 
2131  if ( slp_o<880.e2 .or. slp_o>min(slp_env,slp_mask) .or. abs(pos(2))*rad2deg>40. ) goto 5000 ! next storm
2132 
2133  if ( r_vor < 30.e3 ) then
2134  r_vor = r_min + (slp_env-slp_o)/20.e2*r_inc ! radius of influence
2135  endif
2136 
2137  do j=js, je
2138  do i=is, ie
2139  dist = great_circle_dist(pos, agrid(i,j,1:2), radius)
2140  if( dist < 6.*r_vor ) then
2141  mask(i,j) = mask(i,j) * ( 1. - mask_fac*exp(-(0.5*dist/r_vor)**2)*min(1.,(slp_env-slp_o)/10.e2) )
2142  endif
2143  enddo ! i-loop
2144  enddo ! end j-loop
2145 
2146 5000 continue
2147 
2148  end subroutine get_tc_mask
2149 
2155  subroutine breed_slp_inline(nstep, dt, npz, ak, bk, phis, pe, pk, peln, pkz, delp, u, v, pt, q, nwat, &
2156  zvir, gridstruct, ks, domain_local, bd, hydrostatic)
2157 ! Input
2158  integer, intent(in):: nstep, npz, nwat, ks
2159  real, intent(in):: dt
2160  real, intent(in):: zvir
2161  real, intent(in), dimension(npz+1):: ak, bk
2162  logical, intent(in):: hydrostatic
2163  type(fv_grid_bounds_type), intent(IN) :: bd
2164  real, intent(in):: phis(isd:ied,jsd:jed)
2165  type(domain2d), intent(INOUT) :: domain_local
2166 ! Input/Output
2167  real, intent(inout):: u(isd:ied,jsd:jed+1,npz)
2168  real, intent(inout):: v(isd:ied+1,jsd:jed,npz)
2169  real, intent(inout), dimension(isd:ied,jsd:jed,npz):: delp, pt
2170  real, intent(inout)::q(isd:ied,jsd:jed,npz,*)
2171 
2172  real, intent(inout):: pk(is:ie,js:je, npz+1)
2173  real, intent(inout):: pe(is-1:ie+1, npz+1,js-1:je+1)
2174  real, intent(inout):: pkz(is:ie,js:je,npz)
2175  real, intent(out):: peln(is:ie,npz+1,js:je)
2176 
2177  type(fv_grid_type), target :: gridstruct
2178 ! local
2179  type(time_type):: time
2180  real:: ps(is:ie,js:je)
2181  real:: dist(is:ie,js:je)
2182  real:: tm(is:ie,js:je)
2183  real:: slp(is:ie,js:je)
2184  real(kind=R_GRID):: pos(2)
2185  real:: slp_o ! sea-level pressure (Pa)
2186  real:: w10_o, p_env
2187  real:: r_vor
2188  real:: relx0, relx, f1, pbreed, pbtop, delp0, dp0
2189  real:: ratio, p_count, p_sum, a_sum, mass_sink, delps
2190  real:: p_lo, p_hi, tau_vt, mslp0
2191  real:: split_time, fac, pdep, r2, r3
2192  integer year, month, day, hour, minute, second
2193  integer n, i, j, k, iq, k0
2194 
2195  real, pointer :: area(:,:)
2196  real(kind=R_GRID), pointer :: agrid(:,:,:)
2197 #ifdef MULTI_GASES
2198  real :: kappax(is:ie,js:je,npz)
2199 #endif
2200 
2201  if ( forecast_mode ) return
2202 
2203  agrid => gridstruct%agrid_64
2204  area => gridstruct%area
2205 
2206  if ( nstorms==0 ) then
2207  if(master) write(*,*) 'NO TC data to process'
2208  return
2209  endif
2210 
2211  if ( k_breed==0 ) k_breed = max(1, ks)
2212 
2213 
2214 ! Advance (local) time
2215  call get_date(fv_time, year, month, day, hour, minute, second)
2216  if ( year /= year_track_data ) then
2217  if (master) write(*,*) 'Warning: The year in storm track data is not the same as model year'
2218  return
2219  endif
2220  time = fv_time ! fv_time is the time at past time step (set in fv_diag)
2221  split_time = calday(year, month, day, hour, minute, second) + dt*real(nstep)/86400.
2222 
2223  elapsed_time = elapsed_time + dt
2224  if ( elapsed_time > nudged_time + 0.1 ) then
2225  if(print_end_breed) then
2226  print_end_breed = .false.
2227  if (master) write(*,*) '*** Vortext Breeding Ended at', day, hour, minute, second
2228  endif
2229  return ! time to return to forecast mode
2230  endif
2231 
2232 !$OMP parallel do default(none) shared(is,ie,js,je,npz,ps,ak,delp,tm,pkz,pt)
2233  do j=js,je
2234 ! ---- Compute ps
2235  do i=is,ie
2236  ps(i,j) = ak(1)
2237  enddo
2238  do k=1,npz
2239  do i=is,ie
2240  ps(i,j) = ps(i,j) + delp(i,j,k)
2241  enddo
2242  enddo
2243 ! Compute lowest layer air temperature:
2244  do i=is,ie
2245  tm(i,j) = pkz(i,j,npz)*pt(i,j,npz) / cp_air ! virtual temp
2246  enddo
2247  enddo
2248 ! call prt_maxmin('TM', tm, is, ie, js, je, 0, 1, 1.)
2249 
2250 !$OMP parallel do default(none) shared(k_breed,npz,conserve_mom,is,ie,js,je,u,v,delp, &
2251 !$OMP conserve_hgt,hydrostatic,pt,pkz,q,nwat)
2252  do k=k_breed+1,npz
2253 
2254  if ( conserve_mom ) then
2255  do j=js,je+1
2256  do i=is,ie
2257  u(i,j,k) = u(i,j,k) * (delp(i,j-1,k)+delp(i,j,k))
2258  enddo
2259  enddo
2260  do j=js,je
2261  do i=is,ie+1
2262  v(i,j,k) = v(i,j,k) * (delp(i-1,j,k)+delp(i,j,k))
2263  enddo
2264  enddo
2265  endif
2266 
2267  if ( conserve_hgt .and. hydrostatic ) then
2268  do j=js,je
2269  do i=is,ie
2270 ! Conserve total enthalpy
2271  pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)*delp(i,j,k)
2272  enddo
2273  enddo
2274  endif
2275 
2276 ! Convert tracer moist mixing ratio to mass
2277  do iq=1,nwat
2278  do j=js,je
2279  do i=is,ie
2280  q(i,j,k,iq) = q(i,j,k,iq) * delp(i,j,k)
2281  enddo
2282  enddo
2283  enddo
2284 
2285  enddo
2286 
2287  do 5000 n=1,nstorms ! loop through all storms
2288 
2289  if ( nobs_tc(n) < min_nobs ) goto 5000
2290 
2291 ! Check min MSLP
2292  mslp0 = 1013.e2
2293  do i=1,nobs_tc(n)
2294  mslp0 = min( mslp0, mslp_obs(i,n) )
2295  enddo
2296  if ( mslp0 > min_mslp ) goto 5000
2297 
2298 !----------------------------------------
2299 ! Obtain slp observation
2300 !----------------------------------------
2301  call get_slp_obs(time, nobs_tc(n), x_obs(1,n), y_obs(1,n), wind_obs(1,n), mslp_obs(1,n), mslp_out(1,n), rad_out(1,n), &
2302  time_tc(1,n), pos(1), pos(2), w10_o, slp_o, r_vor, p_env, stime=split_time, fact=fac)
2303 
2304  if ( slp_o<87500. .or. slp_o>slp_env .or. abs(pos(2))*rad2deg>45. ) then
2305  goto 5000 ! next storm
2306  endif
2307 
2308 
2309  if(nudge_debug .and. master) &
2310  write(*,*) 'Vortex breeding for TC:', n, ' slp=',slp_o/100.,pos(1)*rad2deg,pos(2)*rad2deg
2311 
2312 ! Determine pbtop (top pressure of vortex breeding)
2313  k0 = k_breed
2314 
2315  if ( use_high_top ) then
2316 #ifdef HIGH_TEST
2317  if ( slp_o > 1000.e2 ) then
2318  pbtop = 850.e2
2319  else
2320 ! pbtop = max(200.E2, 850.E2-20.*(1000.E2-slp_o))
2321 ! mp87:
2322  pbtop = max(100.e2, 850.e2-25.*(1000.e2-slp_o))
2323 #else
2324  if ( slp_o > 1000.e2 ) then
2325  pbtop = 750.e2
2326  else
2327  pbtop = max(100.e2, 750.e2-20.*(1000.e2-slp_o))
2328 #endif
2329  endif
2330  else
2331 ! Lower top for vrotex breeding
2332  if ( slp_o > 1000.e2 ) then
2333  pbtop = 900.e2
2334  else
2335  pbtop = max(500.e2, 900.e2-5.*(1000.e2-slp_o)) ! mp48
2336  endif
2337  endif
2338 
2339  do k=1,npz
2340  pbreed = ak(k) + bk(k)*1.e5
2341  if ( pbreed>pbtop ) then
2342  k0 = k
2343  exit
2344  endif
2345  enddo
2346  k0 = max(k0, k_breed)
2347 
2348  do j=js, je
2349  do i=is, ie
2350  dist(i,j) = great_circle_dist( pos, agrid(i,j,1:2), radius)
2351  enddo
2352  enddo
2353 
2354  call compute_slp(is, ie, js, je, tm, ps, phis(is:ie,js:je), slp)
2355 
2356  if ( r_vor < 30.e3 .or. p_env<900.e2 ) then
2357 
2358 ! Compute r_vor & p_env
2359  r_vor = r_min + (slp_env-slp_o)/25.e2*r_inc
2360 
2361 123 continue
2362  p_count = 0.
2363  p_sum = 0.
2364  a_sum = 0.
2365  do j=js, je
2366  do i=is, ie
2367  if( dist(i,j)<(r_vor+del_r) .and. dist(i,j)>r_vor .and. phis(i,j)<250.*grav ) then
2368  p_count = p_count + 1.
2369  p_sum = p_sum + slp(i,j)*area(i,j)
2370  a_sum = a_sum + area(i,j)
2371  endif
2372  enddo
2373  enddo
2374 
2375  call mp_reduce_sum(p_count)
2376 
2377  if ( p_count<32. ) then
2378  if(nudge_debug .and. master) write(*,*) p_count, 'Skipping obs: too few p_count'
2379  goto 5000
2380  endif
2381 
2382  call mp_reduce_sum(p_sum)
2383  call mp_reduce_sum(a_sum)
2384  p_env = p_sum / a_sum
2385 
2386  if(nudge_debug .and. master) write(*,*) 'Environmental SLP=', p_env/100., ' computed radius=', r_vor/1.e3
2387 
2388  if ( p_env>1015.e2 .or. p_env<920.e2 ) then
2389  if( nudge_debug ) then
2390  if(master) write(*,*) 'Environmental SLP out of bound; skipping obs. p_count=', p_count, p_sum
2391  call prt_maxmin('SLP_breeding', slp, is, ie, js, je, 0, 1, 0.01)
2392  endif
2393  goto 5000
2394  endif
2395 
2396  endif
2397 
2398 
2399  if ( p_env < max(pre0_env, slp_o + 200.0) ) then
2400 
2401  r_vor = r_vor + 25.e3
2402 
2403  if(nudge_debug .and. master) then
2404  write(*,*) 'Computed environmental SLP too low'
2405  write(*,*) ' ', p_env/100., slp_o/100.,pos(1)*rad2deg, pos(2)*rad2deg
2406  endif
2407 
2408  if ( slp_o > 1003.e2 .and. r_vor > 1000.e3 ) then
2409 ! if(master) write(*,*) 'Failed to size the Vortex for the weak storm'
2410  goto 5000
2411  endif
2412 
2413  if ( r_vor < 1250.e3 ) goto 123
2414 
2415 ! if(master) write(*,*) 'Failed to size the Vortex; skipping this storm'
2416  goto 5000
2417 
2418  endif
2419 
2420  tau_vt = tau_vt_slp * (1. + (960.e2-slp_o)/100.e2 )
2421  tau_vt = max(abs(dt), tau_vt)
2422 
2423  if ( do_adiabatic_init ) then
2424  relx0 = min(1., 2.*abs(dt)/tau_vt)
2425  else
2426  relx0 = min(1., abs(dt)/tau_vt)
2427  endif
2428 
2429  mass_sink = 0.
2430 
2431 !$OMP parallel do default(none) shared(is,ie,js,je,dist,r_vor,phis,p_env,slp_o,r_hi,r_lo, &
2432 !$OMP ps,tm,relx0,tau_vt_rad,dps_min,slp,ak,k0,delp,npz,area) &
2433 !$OMP private(f1, p_hi, p_lo, relx, delps, pbreed, mass_sink, dp0, pdep)
2434  do j=js, je
2435  do i=is, ie
2436  if( dist(i,j) < r_vor .and. phis(i,j)<250.*grav ) then
2437  f1 = dist(i,j)/r_vor
2438 ! Compute p_obs: assuming local radial distributions of slp are Gaussian
2439  p_hi = p_env - (p_env-slp_o) * exp( -r_hi*f1**2 ) ! upper bound
2440  p_lo = p_env - (p_env-slp_o) * exp( -r_lo*f1**2 ) ! lower bound
2441 
2442  if ( ps(i,j) > p_hi .and. tm(i,j) < tm_max ) then
2443 ! do nothing if lowest layer is too hot
2444 ! Under-development:
2445  relx = relx0*exp( -tau_vt_rad*f1**2 )
2446  delps = relx*(ps(i,j) - p_hi) * min( (tm_max-tm(i,j))/10., 1.)
2447 ! After mp115
2448 ! delps = relx*(ps(i,j) - p_hi) * min( (tm_max-tm(i,j))/5., 1.)
2449  ! Cap the increment to prevent overheating
2450  ! Note: ps is used here to prevent
2451  ! over deepening over terrain
2452  delps = min(delps, dps_min)
2453  elseif ( slp(i,j) < p_lo ) then
2454 ! Over-development:
2455  relx = max(0.5, relx0)
2456  delps = relx*(slp(i,j) - p_lo) ! Note: slp is used here
2457  else
2458  goto 400 ! do nothing; proceed to next storm
2459  endif
2460 
2461 #ifdef SIM_TEST
2462  pbreed = ak(1)
2463  do k=1,k0
2464  pbreed = pbreed + delp(i,j,k)
2465  enddo
2466  f1 = 1. - delps/(ps(i,j)-pbreed)
2467  do k=k0+1,npz
2468  delp(i,j,k) = delp(i,j,k)*f1
2469  enddo
2470  mass_sink = mass_sink + delps*area(i,j)
2471 #else
2472  if ( delps > 0. ) then
2473  pbreed = ak(1)
2474  do k=1,k0
2475  pbreed = pbreed + delp(i,j,k)
2476  enddo
2477  f1 = 1. - delps/(ps(i,j)-pbreed)
2478  do k=k0+1,npz
2479  delp(i,j,k) = delp(i,j,k)*f1
2480  enddo
2481  mass_sink = mass_sink + delps*area(i,j)
2482  else
2483  dp0 = abs(delps)
2484  do k=npz,k0+1,-1
2485  if ( abs(delps) < 1. ) then
2486  delp(i,j,k) = delp(i,j,k) - delps
2487  mass_sink = mass_sink + delps*area(i,j)
2488  go to 400
2489  else
2490  pdep = max(1.0, min(abs(0.4*delps), 0.2*dp0, 0.02*delp(i,j,k)))
2491  pdep = - min(pdep, abs(delps))
2492  delp(i,j,k) = delp(i,j,k) - pdep
2493  delps = delps - pdep
2494  mass_sink = mass_sink + pdep*area(i,j)
2495  endif
2496  enddo
2497  endif
2498 #endif
2499 
2500  endif
2501 400 continue
2502  enddo ! end i-loop
2503  enddo ! end j-loop
2504 
2505  call mp_reduce_sum(mass_sink)
2506  if ( abs(mass_sink)<1.e-40 ) goto 5000
2507 
2508  r2 = r_vor + del_r
2509  r3 = min( 4.*r_vor, max(2.*r_vor, 2500.e3) ) + del_r
2510 
2511  p_sum = 0.
2512  do j=js, je
2513  do i=is, ie
2514  if( dist(i,j)<r3 .and. dist(i,j)>r2 ) then
2515  p_sum = p_sum + area(i,j)
2516  endif
2517  enddo
2518  enddo
2519 
2520  call mp_reduce_sum(p_sum)
2521  mass_sink = mass_sink / p_sum ! mean delta pressure to be added back to the environment to conserve mass
2522  if(master .and. nudge_debug) write(*,*) 'TC#',n, 'Mass tele-ported (pa)=', mass_sink
2523 
2524 !$OMP parallel do default(none) shared(is,ie,js,je,dist,r3,r2,ak,k_breed,delp,ps,mass_sink,npz) &
2525 !$OMP private(pbreed, f1)
2526  do j=js, je
2527  do i=is, ie
2528  if( dist(i,j)<r3 .and. dist(i,j)>r2 ) then
2529  pbreed = ak(1)
2530  do k=1,k_breed
2531  pbreed = pbreed + delp(i,j,k)
2532  enddo
2533  f1 = 1. + mass_sink/(ps(i,j)-pbreed)
2534  do k=k_breed+1,npz
2535  delp(i,j,k) = delp(i,j,k)*f1
2536  enddo
2537  endif
2538  enddo
2539  enddo
2540 
2541 ! ---- re-compute ps
2542  do j=js,je
2543  do i=is,ie
2544  ps(i,j) = ak(1)
2545  enddo
2546  do k=1,npz
2547  do i=is,ie
2548  ps(i,j) = ps(i,j) + delp(i,j,k)
2549  enddo
2550  enddo
2551  enddo
2552 
2553 5000 continue
2554 
2555 !--------------------------
2556 ! Update delp halo regions:
2557 !--------------------------
2558  call mpp_update_domains(delp, domain_local, complete=.true.)
2559 
2560 !$OMP parallel do default(none) shared(is,ie,js,je,npz,pe,ak,delp,k_breed,peln,pk)
2561  do j=js-1,je+1
2562  do i=is-1,ie+1
2563  pe(i,1,j) = ak(1)
2564  enddo
2565  do k=2,npz+1
2566  do i=is-1,ie+1
2567  pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
2568  enddo
2569  enddo
2570  enddo
2571 
2572 #ifdef MULTI_GASES
2573 !$OMP parallel do default(none) shared(is,ie,js,je,npz,k_breed,kappax,num_gas)
2574  do k=k_breed,npz
2575  do j=js,je
2576  do i=is,ie
2577  kappax(i,j,k)= virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))
2578  enddo
2579  enddo
2580  enddo
2581 #endif
2582  do k=k_breed+1,npz+1
2583  do j=js,je
2584  do i=is,ie
2585  peln(i,k,j) = log(pe(i,k,j))
2586  pk(i,j,k) = pe(i,k,j)**kappa
2587  enddo
2588  enddo
2589  enddo
2590 
2591 !$OMP parallel do default(none) shared(k_breed,npz,is,ie,js,je,conserve_mom,u,v,delp, &
2592 #ifdef MULTI_GASES
2593 !$OMP kappax, &
2594 #endif
2595 !$OMP conserve_hgt,hydrostatic,pkz,pk,pt,peln)
2596  do k=k_breed+1,npz
2597 
2598  if ( conserve_mom ) then
2599  do j=js,je+1
2600  do i=is,ie
2601  u(i,j,k) = u(i,j,k) / (delp(i,j-1,k)+delp(i,j,k))
2602  enddo
2603  enddo
2604  do j=js,je
2605  do i=is,ie+1
2606  v(i,j,k) = v(i,j,k) / (delp(i-1,j,k)+delp(i,j,k))
2607  enddo
2608  enddo
2609  endif
2610 
2611  if ( conserve_hgt .and. hydrostatic ) then
2612  do j=js,je
2613  do i=is,ie
2614 ! Conserve total enthalpy (static energy)
2615  pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
2616 #ifdef MULTI_GASES
2617  pkz(i,j,k) = exp(kappax(i,j,k)*log(pkz(i,j,k)))
2618 #endif
2619  pt(i,j,k) = pt(i,j,k) / (pkz(i,j,k)*delp(i,j,k))
2620  enddo
2621  enddo
2622 ! pkz last used
2623  endif
2624  enddo
2625 
2626 ! Convert tracer mass back to moist mixing ratio
2627 
2628 !$OMP parallel do default(none) shared(nwat,k_breed,npz,is,ie,js,je,q,delp)
2629  do iq=1,nwat
2630  do k=k_breed+1,npz
2631  do j=js,je
2632  do i=is,ie
2633  q(i,j,k,iq) = q(i,j,k,iq) / delp(i,j,k)
2634  enddo
2635  enddo
2636  enddo
2637  enddo
2638 
2639  if(hydrostatic) call mpp_update_domains(pt, domain_local, complete=.true.)
2640 
2641  nullify(agrid)
2642  nullify(area)
2643 
2644  end subroutine breed_slp_inline
2645 
2648  subroutine breed_srf_w10(time, dt, npz, ak, bk, ps, phis, slp, delp, u, v, gridstruct)
2649  type(time_type), intent(in):: time
2650  integer, intent(in):: npz
2651  real, intent(in):: dt
2652  real, intent(in), dimension(npz+1):: ak, bk
2653  real, intent(in):: phis(isd:ied,jsd:jed)
2654  real, intent(in):: ps(isd:ied,jsd:jed)
2655  real, intent(in), dimension(is:ie,js:je):: slp
2656  type(fv_grid_type), intent(IN), target :: gridstruct
2657  real, intent(inout):: u(isd:ied,jsd:jed+1,npz)
2658  real, intent(inout):: v(isd:ied+1,jsd:jed,npz)
2659 ! Input/Output
2660  real, intent(inout), dimension(isd:ied,jsd:jed,npz):: delp
2661 ! local
2662  real, dimension(is:ie,js:je):: us, vs
2663  real wu(is:ie, js:je+1)
2664  real wv(is:ie+1,js:je)
2665  real u1(is:ie), v1(is:ie)
2666  real:: dist(is:ie,js:je), wind(is:ie,js:je)
2667  real(kind=R_GRID):: pos(2)
2668  real:: slp_o
2669  real:: w10_o, p_env
2670  real:: r_vor, pc, p_count
2671  real:: r_max, speed, ut, vt, speed_local ! tangent wind speed
2672  real:: u_bg, v_bg, mass, t_mass
2673  real:: relx0, relx, f1
2674  real:: z0, mslp0
2675  real:: zz = 35.
2676  real:: wind_fac
2677  integer n, i, j
2678 
2679  ! Pointers
2680  real, pointer, dimension(:,:) :: dx, dy, rdxa, rdya, a11, a21, a12, a22, area
2681  real(kind=R_GRID), pointer, dimension(:,:,:) :: agrid, vlon, vlat
2682 
2683  dx => gridstruct%dx
2684  dy => gridstruct%dy
2685  rdxa => gridstruct%rdxa
2686  rdya => gridstruct%rdya
2687  a11 => gridstruct%a11
2688  a21 => gridstruct%a21
2689  a12 => gridstruct%a12
2690  a22 => gridstruct%a22
2691  area => gridstruct%area
2692  agrid => gridstruct%agrid_64
2693  vlon => gridstruct%vlon
2694  vlat => gridstruct%vlat
2695 
2696 
2697  if ( nstorms==0 ) then
2698  if(master) write(*,*) 'NO TC data to process'
2699  return
2700  endif
2701 
2702 ! Compute lat-lon winds on A grid
2703 !$OMP parallel do default(none) shared(is,ie,js,je,wu,u,npz,dx)
2704  do j=js,je+1
2705  do i=is,ie
2706  wu(i,j) = u(i,j,npz)*dx(i,j)
2707  enddo
2708  enddo
2709 !$OMP parallel do default(none) shared(is,ie,js,je,wv,v,npz,dy)
2710  do j=js,je
2711  do i=is,ie+1
2712  wv(i,j) = v(i,j,npz)*dy(i,j)
2713  enddo
2714  enddo
2715 
2716 !$OMP parallel do default(none) shared(is,ie,js,je,wu,rdxa,wv,rdya,us,vs,a11,a12,a21,a22,wind) &
2717 !$OMP private(u1,v1)
2718  do j=js, je
2719  do i=is, ie
2720 ! Co-variant to Co-variant "vorticity-conserving" interpolation
2721  u1(i) = (wu(i,j) + wu(i,j+1)) * rdxa(i,j)
2722  v1(i) = (wv(i,j) + wv(i+1,j)) * rdya(i,j)
2723 ! Cubed (cell center co-variant winds) to lat-lon:
2724  us(i,j) = a11(i,j)*u1(i) + a12(i,j)*v1(i)
2725  vs(i,j) = a21(i,j)*u1(i) + a22(i,j)*v1(i)
2726 ! Wind speed
2727  wind(i,j) = sqrt( us(i,j)**2 + vs(i,j)**2 )
2728  enddo
2729  enddo
2730 
2731  relx0 = min(1., dt/tau_vt_wind)
2732 
2733  do 3000 n=1,nstorms ! loop through all storms
2734 
2735  if ( nobs_tc(n) < min_nobs ) goto 3000
2736 
2737 ! Check min MSLP
2738  mslp0 = 1013.e2
2739  do i=1,nobs_tc(n)
2740  mslp0 = min( mslp0, mslp_obs(i,n) )
2741  enddo
2742  if ( mslp0 > min_mslp ) goto 3000
2743 
2744 !----------------------------------------
2745 ! Obtain slp observation
2746 !----------------------------------------
2747  call get_slp_obs(time, nobs_tc(n), x_obs(1,n), y_obs(1,n), wind_obs(1,n), mslp_obs(1,n), mslp_out(1,n), rad_out(1,n), &
2748  time_tc(1,n), pos(1), pos(2), w10_o, slp_o, r_vor, p_env)
2749 
2750  if ( slp_o<90000. .or. slp_o>slp_env .or. abs(pos(2))*rad2deg>35. ) goto 3000 ! next storm
2751 
2752 
2753  do j=js, je
2754  do i=is, ie
2755  dist(i,j) = great_circle_dist( pos, agrid(i,j,1:2), radius )
2756  enddo
2757  enddo
2758 
2759  r_vor = r_min + (slp_env-slp_o)/25.e2*r_inc
2760 
2761 !----------------------------------------------------
2762 ! * Find model's SLP center nearest to the observation
2763 ! * Find maximum wind speed at the lowest model level
2764 
2765  speed_local = 0.
2766  r_max = -999.
2767  pc = 1013.e2
2768  p_count = 0.
2769  do j=js, je
2770  do i=is, ie
2771  if( dist(i,j) < r_vor ) then
2772 
2773 ! Counting the "land" points:
2774  if ( dist(i,j)<0.5*r_vor .and. abs(phis(i,j))>2.*grav ) p_count = p_count + 1.
2775 
2776  pc = min(pc, slp(i,j))
2777 
2778  if ( speed_local < wind(i,j) ) then
2779  speed_local = wind(i,j)
2780  r_max = dist(i,j)
2781  endif
2782 
2783  endif
2784  enddo
2785  enddo
2786 
2787  call mp_reduce_sum(p_count)
2788  if ( p_count>32 ) goto 3000 ! over/near rough land
2789 
2790  if ( w10_o < 0. ) then ! 10-m wind obs is not available
2791 ! Uses Atkinson_Holliday wind-pressure correlation
2792  w10_o = 3.446778 * (1010.-slp_o/100.)**0.644
2793  endif
2794 
2795  speed = speed_local
2796  call mp_reduce_max(speed) ! global max wind (near storm)
2797  call mp_reduce_min(pc)
2798 
2799  if ( speed_local < speed ) then
2800  r_max = -999.
2801  endif
2802  call mp_reduce_max(r_max)
2803  if( r_max<0. ) call mpp_error(fatal,'==> Error in r_max')
2804 
2805 ! ---------------------------------------------------
2806 ! Determine surface wind speed and radius for nudging
2807 ! ---------------------------------------------------
2808 
2809 ! Compute surface roughness z0 from w10, based on Eq (4) & (5) from Moon et al. 2007
2810  if ( w10_o > 12.5 ) then
2811  z0 = (0.085*w10_o - 0.58) * 1.e-3
2812 ! z0 (w10=40) = 2.82E-3
2813 ! z0 = min( z0, 2.82E-3 )
2814 
2815  else
2816  z0 = 0.0185/grav*(0.001*w10_o**2 + 0.028*w10_o)**2
2817  endif
2818 
2819 ! lowest layer height: zz
2820 
2821  wind_fac = log(zz/z0) / log(10./z0)
2822  if( nudge_debug .and. master ) write(*,*) 'Wind adjustment factor=', wind_fac
2823  if( wind_fac<1. ) call mpp_error(fatal,'==> Error in wind_fac')
2824 
2825  if ( pc < slp_o ) then
2826 !--
2827 ! The storm in the model is over developed
2828 ! if ( (pc+3.0E2)>slp_o .or. speed <= w10_o ) go to 3000 ! next storm
2829 !--
2830 ! using radius (r_max) as dtermined above;
2831 ! What if model's pressure center is very far from the observed?
2832 ! using obs wind
2833  speed = wind_fac*w10_o
2834  else
2835 ! The storm in the model is under developed; using max wind
2836  speed = max(wind_fac*w10_o, speed)
2837  if ( pc>1009.e2 ) r_max = 0.5 * r_vor
2838  endif
2839 
2840 ! More adjustment
2841 
2842 ! Some bounds on the radius of maximum wind:
2843  r_max = max(2.5*grid_size, r_max) ! at least 2.5X the grid size
2844  r_max = min(0.75*r_vor, r_max)
2845 
2846  t_mass = 0.
2847  u_bg = 0.
2848  v_bg = 0.
2849 
2850  if ( add_bg_wind ) then
2851  p_count = 0.
2852  do j=js, je
2853  do i=is, ie
2854  if( dist(i,j) <= min(r_vor,2.*r_fac*r_max) .and. phis(i,j)<2.*grav ) then
2855  mass = area(i,j)*delp(i,j,npz)
2856 !-- using model winds ----------------------------------
2857  u_bg = u_bg + us(i,j)*mass
2858  v_bg = v_bg + vs(i,j)*mass
2859 !-------------------------------------------------------
2860  t_mass = t_mass + mass
2861  p_count = p_count + 1.
2862  endif
2863  enddo
2864  enddo
2865  call mp_reduce_sum(p_count)
2866  if ( p_count<16. ) go to 3000
2867 
2868  call mp_reduce_sum(t_mass)
2869  call mp_reduce_sum(u_bg)
2870  call mp_reduce_sum(v_bg)
2871  u_bg = u_bg / t_mass
2872  v_bg = v_bg / t_mass
2873  if ( nudge_debug .and. master ) write(*,*) pos(2)*rad2deg, 'vortex bg wind=', u_bg, v_bg
2874  endif
2875 
2876  relx = relx0
2877 ! Nudge wind in the "inner core":
2878 #ifdef TTT
2879  do j=js, je
2880  do i=is, ie
2881  if( dist(i,j) <= min(r_vor, r_fac*r_max) .and. phis(i,j)<2.*grav ) then
2882  f1 = dist(i,j)/r_max
2883  relx = relx0*exp( -tau_vt_rad*f1**2 )
2884  if( dist(i,j)<=r_max ) then
2885  speed_local = speed * f1
2886  else
2887  speed_local = speed / f1**0.75
2888  endif
2889  call tangent_wind(vlon(i,j,1:3), vlat(i,j,1:3), speed_local, pos, agrid(i,j,1:2), ut, vt)
2890  ut = ut + u_bg
2891  vt = vt + v_bg
2892 ! Update:
2893  us(i,j) = relx*(ut-us(i,j))
2894  vs(i,j) = relx*(vt-vs(i,j))
2895  endif
2896 400 continue
2897  enddo ! end i-loop
2898  enddo ! end j-loop
2899 #else
2900  do j=js, je
2901  do i=is, ie
2902  if( dist(i,j) <= min(r_vor, r_fac*r_max) .and. phis(i,j)<2.*grav ) then
2903  f1 = dist(i,j)/r_max
2904  relx = relx0*exp( -tau_vt_rad*f1**2 )
2905  if( dist(i,j)<=r_max ) then
2906  speed_local = speed * f1
2907  else
2908  speed_local = speed / f1**0.75
2909  endif
2910  call tangent_wind(vlon(i,j,1:3), vlat(i,j,1:3), speed_local, pos, agrid(i,j,1:2), ut, vt)
2911  ut = ut + u_bg
2912  vt = vt + v_bg
2913 ! Update:
2914  us(i,j) = relx*(ut-us(i,j))
2915  vs(i,j) = relx*(vt-vs(i,j))
2916  endif
2917 400 continue
2918  enddo ! end i-loop
2919  enddo ! end j-loop
2920 #endif
2921 
2922 3000 continue
2923 
2924  end subroutine breed_srf_w10
2925 
2927  subroutine breed_srf_winds(time, dt, npz, u_obs, v_obs, ak, bk, ps, phis, delp, ua, va, u_dt, v_dt, pt, q, nwat, zvir, gridstruct)
2928 ! Input
2929  type(time_type), intent(in):: time
2930  integer, intent(in):: npz, nwat
2931  real, intent(in):: dt
2932  real, intent(in):: zvir
2933  real, intent(in), dimension(npz+1):: ak, bk
2934  real, intent(in), dimension(isd:ied,jsd:jed):: phis, ps
2935  real, intent(in), dimension(is:ie,js:je,npz):: u_obs, v_obs
2936  type(fv_grid_type), intent(IN), target :: gridstruct
2937 ! Input/Output
2938  real, intent(inout), dimension(isd:ied,jsd:jed,npz):: delp, pt, ua, va, u_dt, v_dt
2939  real, intent(inout)::q(isd:ied,jsd:jed,npz,nwat)
2940 ! local
2941  real:: dist(is:ie,js:je), wind(is:ie,js:je)
2942  real:: slp(is:ie,js:je)
2943  real(kind=R_GRID):: pos(2)
2944  real:: slp_o
2945  real:: w10_o, p_env
2946  real:: r_vor, pc, p_count
2947  real:: r_max, speed, ut, vt, speed_local ! tangent wind speed
2948  real:: u_bg, v_bg, mass, t_mass
2949  real:: relx0, relx, f1, rdt
2950  real:: z0, mslp0
2951  real:: zz = 35.
2952  real:: wind_fac
2953  integer n, i, j, k, iq
2954 
2955  real, pointer :: area(:,:)
2956  real(kind=R_GRID), pointer :: vlon(:,:,:), vlat(:,:,:), agrid(:,:,:)
2957 
2958  area => gridstruct%area
2959  vlon => gridstruct%vlon
2960  vlat => gridstruct%vlat
2961  agrid => gridstruct%agrid_64
2962 
2963  if ( nstorms==0 ) then
2964  if(master) write(*,*) 'NO TC data to process'
2965  return
2966  endif
2967 
2968  rdt = 1./dt
2969  relx0 = min(1., dt/tau_vt_wind)
2970 
2971 !$OMP parallel do default(none) shared(is,ie,js,je,slp,ps,phis,pt,npz,wind,ua,va)
2972  do j=js, je
2973  do i=is, ie
2974  slp(i,j) = ps(i,j)*exp(phis(i,j)/(rdgas*(pt(i,j,npz)+3.25e-3*phis(i,j)/grav)))
2975  wind(i,j) = sqrt( ua(i,j,npz)**2 + va(i,j,npz)**2 )
2976  enddo
2977  enddo
2978 
2979 !!!!!$OMP parallel do default(none) private(pos, w10_o, slp_o, r_vor, p_env)
2980  do 3000 n=1,nstorms ! loop through all storms
2981 
2982  if ( nobs_tc(n) < min_nobs ) goto 3000
2983 
2984 ! Check min MSLP
2985  mslp0 = 1013.e2
2986  do i=1,nobs_tc(n)
2987  mslp0 = min( mslp0, mslp_obs(i,n) )
2988  enddo
2989  if ( mslp0 > min_mslp ) goto 3000
2990 
2991 !----------------------------------------
2992 ! Obtain slp observation
2993 !----------------------------------------
2994  call get_slp_obs(time, nobs_tc(n), x_obs(1,n), y_obs(1,n), wind_obs(1,n), mslp_obs(1,n), mslp_out(1,n), rad_out(1,n), &
2995  time_tc(1,n), pos(1), pos(2), w10_o, slp_o, r_vor, p_env)
2996 
2997  if ( slp_o<90000. .or. slp_o>slp_env .or. abs(pos(2))*rad2deg>35. ) goto 3000 ! next storm
2998 
2999 
3000  do j=js, je
3001  do i=is, ie
3002  dist(i,j) = great_circle_dist( pos, agrid(i,j,1:2), radius )
3003  enddo
3004  enddo
3005 
3006  r_vor = r_min + (slp_env-slp_o)/25.e2*r_inc
3007 
3008 !----------------------------------------------------
3009 
3010 ! * Find model's SLP center nearest to the observation
3011 ! * Find maximum wind speed at the lowest model level
3012 
3013  speed_local = 0.
3014  r_max = -999.
3015  pc = 1013.e2
3016  p_count = 0.
3017  do j=js, je
3018  do i=is, ie
3019  if( dist(i,j) < r_vor ) then
3020 
3021 ! Counting the "land" points:
3022  if ( dist(i,j)<0.5*r_vor .and. abs(phis(i,j))>2.*grav ) p_count = p_count + 1.
3023 
3024  pc = min(pc, slp(i,j))
3025 
3026  if ( speed_local < wind(i,j) ) then
3027  speed_local = wind(i,j)
3028  r_max = dist(i,j)
3029  endif
3030 
3031  endif
3032  enddo
3033  enddo
3034 
3035  call mp_reduce_sum(p_count)
3036  if ( p_count>32 ) goto 3000 ! over/near rough land
3037 
3038  if ( w10_o < 0. ) then ! 10-m wind obs is not available
3039 ! Uses Atkinson_Holliday wind-pressure correlation
3040  w10_o = 3.446778 * (1010.-slp_o/100.)**0.644
3041  endif
3042 
3043  speed = speed_local
3044  call mp_reduce_max(speed) ! global max wind (near storm)
3045  call mp_reduce_min(pc)
3046 
3047  if ( speed_local < speed ) then
3048  r_max = -999.
3049  endif
3050  call mp_reduce_max(r_max)
3051  if( r_max<0. ) call mpp_error(fatal,'==> Error in r_max')
3052 
3053 ! ---------------------------------------------------
3054 ! Determine surface wind speed and radius for nudging
3055 ! ---------------------------------------------------
3056 
3057 ! Compute surface roughness z0 from w10, based on Eq (4) & (5) from Moon et al. 2007
3058  if ( w10_o > 12.5 ) then
3059  z0 = (0.085*w10_o - 0.58) * 1.e-3
3060 ! z0 (w10=40) = 2.82E-3
3061 ! z0 = min( z0, 2.82E-3 )
3062 
3063  else
3064  z0 = 0.0185/grav*(0.001*w10_o**2 + 0.028*w10_o)**2
3065  endif
3066 
3067 ! lowest layer height: zz
3068 
3069  wind_fac = log(zz/z0) / log(10./z0)
3070  if( nudge_debug .and. master ) write(*,*) 'Wind adjustment factor=', wind_fac
3071  if( wind_fac<1. ) call mpp_error(fatal,'==> Error in wind_fac')
3072 
3073  if ( pc < slp_o ) then
3074 !--
3075 ! The storm in the model is over developed
3076 ! if ( (pc+3.0E2)>slp_o .or. speed <= w10_o ) go to 3000 ! next storm
3077 !--
3078 ! using radius (r_max) as dtermined above;
3079 ! What if model's pressure center is very far from the observed?
3080 ! using obs wind
3081  speed = wind_fac*w10_o
3082  else
3083 ! The storm in the model is under developed; using max wind
3084  speed = max(wind_fac*w10_o, speed)
3085  if ( pc>1009.e2 ) r_max = 0.5 * r_vor
3086  endif
3087 
3088 ! More adjustment
3089 
3090 ! Some bounds on the radius of maximum wind:
3091  r_max = max(2.5*grid_size, r_max) ! at least 2.5X the grid size
3092  r_max = min(0.75*r_vor, r_max)
3093 
3094  t_mass = 0.
3095  u_bg = 0.
3096  v_bg = 0.
3097 
3098  if ( add_bg_wind ) then
3099  p_count = 0.
3100  do j=js, je
3101  do i=is, ie
3102  if( dist(i,j) <= min(r_vor,2.*r_fac*r_max) .and. phis(i,j)<2.*grav ) then
3103  mass = area(i,j)*delp(i,j,npz)
3104 !-- using model winds ----------------------------------
3105 ! u_bg = u_bg + ua(i,j,npz)*mass
3106 ! v_bg = v_bg + va(i,j,npz)*mass
3107 !-------------------------------------------------------
3108 ! Using analysis winds
3109  u_bg = u_bg + u_obs(i,j,npz)*mass
3110  v_bg = v_bg + v_obs(i,j,npz)*mass
3111  t_mass = t_mass + mass
3112  p_count = p_count + 1.
3113  endif
3114  enddo
3115  enddo
3116  call mp_reduce_sum(p_count)
3117  if ( p_count<16. ) go to 3000
3118 
3119  call mp_reduce_sum(t_mass)
3120  call mp_reduce_sum(u_bg)
3121  call mp_reduce_sum(v_bg)
3122  u_bg = u_bg / t_mass
3123  v_bg = v_bg / t_mass
3124 ! if ( master ) write(*,*) pos(2)*rad2deg, 'vortex bg wind=', u_bg, v_bg
3125  endif
3126 
3127  relx = relx0
3128  k = npz ! lowest layer only
3129 ! Nudge wind in the "inner core":
3130  do j=js, je
3131  do i=is, ie
3132  if( dist(i,j) <= min(r_vor, r_fac*r_max) .and. phis(i,j)<2.*grav ) then
3133  f1 = dist(i,j)/r_max
3134  relx = relx0*exp( -tau_vt_rad*f1**2 )
3135  if( dist(i,j)<=r_max ) then
3136  speed_local = speed * f1
3137  else
3138  speed_local = speed / f1**0.75
3139  endif
3140  call tangent_wind(vlon(i,j,1:3), vlat(i,j,1:3), speed_local, pos, agrid(i,j,1:2), ut, vt)
3141  ut = ut + u_bg
3142  vt = vt + v_bg
3143  u_dt(i,j,k) = u_dt(i,j,k) + relx*(ut-ua(i,j,k)) * rdt
3144  v_dt(i,j,k) = v_dt(i,j,k) + relx*(vt-va(i,j,k)) * rdt
3145 ! Update:
3146  ua(i,j,k) = ua(i,j,k) + relx*(ut-ua(i,j,k))
3147  va(i,j,k) = va(i,j,k) + relx*(vt-va(i,j,k))
3148  endif
3149 400 continue
3150  enddo ! end i-loop
3151  enddo ! end j-loop
3152 
3153 3000 continue
3154 
3155  end subroutine breed_srf_winds
3156 
3157  subroutine tangent_wind ( elon, elat, speed, po, pp, ut, vt )
3158  real, intent(in):: speed
3159  real(kind=R_GRID), intent(in):: po(2), pp(2)
3160  real(kind=R_GRID), intent(in):: elon(3), elat(3)
3161  real, intent(out):: ut, vt
3162 ! local
3163  real(kind=R_GRID):: e1(3), eo(3), ep(3), op(3)
3164 
3165  call latlon2xyz(po, eo)
3166  call latlon2xyz(pp, ep)
3167 
3168  op(:) = ep(:) - eo(:)
3169  call normalize_vect( op )
3170 
3171  call vect_cross(e1, ep, eo)
3172 
3173  ut = speed * (e1(1)*elon(1) + e1(2)*elon(2) + e1(3)*elon(3))
3174  vt = speed * (e1(1)*elat(1) + e1(2)*elat(2) + e1(3)*elat(3))
3175 
3176 ! SH:
3177  if ( po(2) < 0. ) then
3178  ut = -ut
3179  vt = -vt
3180  endif
3181 
3182  end subroutine tangent_wind
3183 
3184 
3185  subroutine get_slp_obs(time, nobs, lon_obs, lat_obs, w10, mslp, slp_out, r_out, time_obs, &
3186  x_o, y_o, w10_o, slp_o, r_vor, p_vor, stime, fact)
3187 ! Input
3188  type(time_type), intent(in):: time
3189  integer, intent(in):: nobs
3190  real(KIND=4), intent(in):: lon_obs(nobs)
3191  real(KIND=4), intent(in):: lat_obs(nobs)
3192  real(KIND=4), intent(in):: w10(nobs)
3193  real(KIND=4), intent(in):: mslp(nobs)
3194  real(KIND=4), intent(in):: slp_out(nobs)
3195  real(KIND=4), intent(in):: r_out(nobs)
3196  real(KIND=4), intent(in):: time_obs(nobs)
3197  real, optional, intent(in):: stime
3198  real, optional, intent(out):: fact
3199 ! Output
3200  real(kind=R_GRID), intent(out):: x_o , y_o
3201  real, intent(out):: w10_o
3202  real, intent(out):: slp_o
3203  real, intent(out):: r_vor, p_vor
3204 ! Internal:
3205  real:: t_thresh
3206  real(kind=R_GRID):: p1(2), p2(2)
3207  real time_model
3208  real(kind=R_GRID) fac
3209  integer year, month, day, hour, minute, second, n
3210 
3211  t_thresh = 600./86400. ! unit: days
3212 
3213  w10_o = -100000.
3214  slp_o = -100000.
3215  x_o = -100.*pi
3216  y_o = -100.*pi
3217  p_vor = -1.e10
3218  r_vor = -1.e10
3219 
3220  if ( present(stime) ) then
3221  time_model = stime
3222  else
3223  call get_date(time, year, month, day, hour, minute, second)
3224 
3225  if ( year /= year_track_data ) then
3226  if (master) write(*,*) 'Warning: The year in storm track data is not the same as model year'
3227  return
3228  endif
3229 
3230  time_model = calday(year, month, day, hour, minute, second)
3231 ! if(nudge_debug .and. master) write(*,*) 'Model:', time_model, year, month, day, hour, minute, second
3232  endif
3233 
3234 !-------------------------------------------------------------------------------------------
3235 ! if ( time_model <= time_obs(1) .or. time_model >= time_obs(nobs) ) then
3236 ! return
3237 !-------------------------------------------------------------------------------------------
3238 
3239  if ( time_model <= (time_obs(1)-t_thresh) .or. time_model >= time_obs(nobs) ) return
3240 
3241  if ( time_model <= time_obs(1) ) then
3242 !--
3243 ! This is an attempt to perform vortex breeding several minutes before the first available observation
3244 !--
3245  w10_o = w10(1)
3246  slp_o = mslp(1)
3247  x_o = lon_obs(1)
3248  y_o = lat_obs(1)
3249  if ( present(fact) ) fact = 1.25
3250  else
3251  do n=1,nobs-1
3252  if( time_model >= time_obs(n) .and. time_model <= time_obs(n+1) ) then
3253  fac = (time_model-time_obs(n)) / (time_obs(n+1)-time_obs(n))
3254  w10_o = w10(n) + ( w10(n+1)- w10(n)) * fac
3255  slp_o = mslp(n) + ( mslp(n+1)- mslp(n)) * fac
3256 ! Trajectory interpolation:
3257 ! Linear in (lon,lat) space
3258 ! x_o = lon_obs(n) + (lon_obs(n+1)-lon_obs(n)) * fac
3259 ! y_o = lat_obs(n) + (lat_obs(n+1)-lat_obs(n)) * fac
3260  p1(1) = lon_obs(n); p1(2) = lat_obs(n)
3261  p2(1) = lon_obs(n+1); p2(2) = lat_obs(n+1)
3262  call intp_great_circle(fac, p1, p2, x_o, y_o)
3263 !----------------------------------------------------------------------
3264  if ( present(fact) ) fact = 1. + 0.25*cos(fac*2.*pi)
3265 ! Additional data from the extended best track
3266 ! if ( slp_out(n)>0. .and. slp_out(n+1)>0. .and. r_out(n)>0. .and. r_out(n+1)>0. ) then
3267 ! p_vor = slp_out(n) + ( slp_out(n+1) - slp_out(n)) * fac
3268 ! r_vor = r_out(n) + ( r_out(n+1) - r_out(n)) * fac
3269 ! endif
3270  return
3271  endif
3272  enddo
3273  endif
3274 
3275  end subroutine get_slp_obs
3276 
3277 
3278  subroutine slp_obs_init
3279  integer:: unit, n, nobs
3280  character(len=3):: GMT
3281  character(len=9):: ts_name
3282  character(len=19):: comment
3283  integer:: mmddhh, yr, year, month, day, hour, MPH, islp
3284  integer:: it, i1, i2, p_ring, d_ring
3285  real:: lon_deg, lat_deg, cald, slp, mps
3286 
3287  nobs_tc(:) = 0
3288  time_tc(:,:) = 0.
3289  wind_obs(:,:) = -100000.
3290  mslp_obs(:,:) = -100000.
3291  x_obs(:,:) = - 100.*pi
3292  y_obs(:,:) = - 100.*pi
3293 
3294  mslp_out(:,:) = -1.e10
3295  rad_out(:,:) = -1.e10
3296 
3297  if( track_file_name == "No_File_specified" ) then
3298  if(master) write(*,*) 'No TC track file specified'
3299  return
3300  else
3301  unit = get_unit()
3302  open( unit, file=track_file_name)
3303  endif
3304 
3305  read(unit, *) year
3306  if(master) write(*,*) 'Reading TC track data for YEAR=', year
3307 
3308  year_track_data = year
3309 
3310  nstorms = 0
3311  nobs = 0
3312  month = 99
3313 
3314  if ( ibtrack ) then
3315 
3316 !---------------------------------------------------------------
3317 ! The data format is from Ming Zhoa's processed ibTrack datasets
3318 !---------------------------------------------------------------
3319 
3320  read(unit, *) ts_name, nobs, yr, month, day, hour
3321 
3322  if ( yr /= year ) then
3323  if(master) write(*, *) 'Year inconsistency found !!!'
3324  call mpp_error(fatal,'==> Error in reading best track data')
3325  endif
3326 
3327  do while ( ts_name=='start' )
3328 
3329  nstorms = nstorms + 1
3330  nobs_tc(nstorms) = nobs ! observation count for this storm
3331  if(nudge_debug .and. master) write(*, *) 'Read Data for TC#', nstorms, nobs
3332 
3333  do it=1, nobs
3334  read(unit, *) lon_deg, lat_deg, mps, slp, yr, month, day, hour
3335 ! if ( yr /= year ) then
3336 ! if(master) write(*, *) 'Extended to year + 1', yr
3337 ! endif
3338  cald = calday(yr, month, day, hour, 0, 0)
3339  time_tc(it,nstorms) = cald
3340  if(nudge_debug .and. master) write(*, 100) cald, month, day, hour, lon_deg, lat_deg, mps, slp
3341 
3342  wind_obs(it,nstorms) = mps ! m/s
3343  mslp_obs(it,nstorms) = 100.*slp
3344  y_obs(it,nstorms) = lat_deg * deg2rad
3345  x_obs(it,nstorms) = lon_deg * deg2rad
3346  enddo
3347 
3348  read(unit, *) ts_name, nobs, yr, month, day, hour
3349  enddo
3350 100 format(1x, f9.2, 1x, i3, 1x, i2, 1x, i2, 1x, f6.1, 1x, f6.1, 1x, f4.1, 1x, f6.1)
3351 
3352  else
3353 
3354  do while ( month /= 0 )
3355 
3356  read(unit, *) month, day, hour, gmt, lat_deg, lon_deg, mph, islp, comment
3357 
3358  select case (month)
3359 
3360  case (99) ! New storm record to start
3361  nstorms = nstorms + 1
3362  nobs = 0
3363  if(master) write(*, *) 'Reading data for TC#', nstorms, comment
3364  case ( 0) ! end of record
3365  if(master) write(*, *) 'End of record reached'
3366  case default
3367  nobs = nobs + 1
3368  cald = calday(year, month, day, hour, 0, 0)
3369  time_tc(nobs,nstorms) = cald
3370  nobs_tc(nstorms) = nobs ! observation count for this storm
3371 
3372  if(master) write(*, 200) nobs, cald, month, day, hour, gmt, lat_deg, lon_deg, mph, islp, comment
3373  mslp_obs(nobs,nstorms) = 100. * real(islp)
3374  y_obs(nobs,nstorms) = lat_deg * deg2rad
3375  if ( gmt == 'GMT' ) then
3376 ! Transfrom x from (-180 , 180) to (0, 360) then to radian
3377  if ( lon_deg < 0 ) then
3378  x_obs(nobs,nstorms) = (360.+lon_deg) * deg2rad
3379  else
3380  x_obs(nobs,nstorms) = (360.-lon_deg) * deg2rad
3381  endif
3382  elseif ( gmt == 'PAC' ) then ! Pacific storms
3383  x_obs(nobs,nstorms) = lon_deg * deg2rad
3384  endif
3385  end select
3386 
3387  enddo
3388 
3389  endif
3390 
3391  close(unit)
3392 
3393  if(master) then
3394  write(*,*) 'TC vortex breeding: total storms=', nstorms
3395  if ( nstorms/=0 ) then
3396  do n=1,nstorms
3397  write(*, *) 'TC#',n, ' contains ', nobs_tc(n),' observations'
3398  enddo
3399  endif
3400  endif
3401 
3402 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)
3403 
3404  end subroutine slp_obs_init
3405 
3407  ! Julian day (0 to 365 for non-leap year)
3408  real function calday(year, month, day, hour, minute, sec)
3409 ! input:
3410  integer, intent(in):: year, month, day, hour
3411  integer, intent(in):: minute, sec
3412 ! Local:
3413  integer n, m, ds, nday
3414  real tsec
3415  integer days(12)
3416  data days /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
3417 
3418  ds = day - 1
3419 
3420  if( month /= 1 ) then
3421  do m=1, month-1
3422  if( m==2 .and. leap_year(year) ) then
3423  ds = ds + 29
3424  else
3425  ds = ds + days(m)
3426  endif
3427  enddo
3428  endif
3429 
3430  if ( leap_year(year_track_data) ) then
3431  nday = 366
3432  else
3433  nday = 365
3434  endif
3435 
3436  calday = real((year-year_track_data)*nday + ds) + real(hour*3600 + minute*60 + sec)/86400.
3437 
3438  end function calday
3439 
3442  logical function leap_year(ny)
3443  integer, intent(in):: ny
3444  integer ny00
3445 !
3446 ! No leap years prior to 0000
3447 !
3448  parameter( ny00 = 0000 )
3449 
3450  if( ny >= ny00 ) then
3451  if( mod(ny,100) == 0. .and. mod(ny,400) == 0. ) then
3452  leap_year = .true.
3453  elseif( mod(ny,4) == 0. .and. mod(ny,100) /= 0. ) then
3454  leap_year = .true.
3455  else
3456  leap_year = .false.
3457  endif
3458  else
3459  leap_year = .false.
3460  endif
3461 
3462  end function leap_year
3463 
3464 
3465  subroutine pmaxmin( qname, a, imax, jmax, fac )
3467  character(len=*) qname
3468  integer imax, jmax
3469  integer i, j
3470  real a(imax,jmax)
3471 
3472  real qmin(jmax), qmax(jmax)
3473  real pmax, pmin
3474  real fac ! multiplication factor
3475 
3476  do j=1,jmax
3477  pmax = a(1,j)
3478  pmin = a(1,j)
3479  do i=2,imax
3480  pmax = max(pmax, a(i,j))
3481  pmin = min(pmin, a(i,j))
3482  enddo
3483  qmax(j) = pmax
3484  qmin(j) = pmin
3485  enddo
3486 !
3487 ! Now find max/min of amax/amin
3488 !
3489  pmax = qmax(1)
3490  pmin = qmin(1)
3491  do j=2,jmax
3492  pmax = max(pmax, qmax(j))
3493  pmin = min(pmin, qmin(j))
3494  enddo
3495 
3496  write(*,*) qname, ' max = ', pmax*fac, ' min = ', pmin*fac
3497 
3498  end subroutine pmaxmin
3499 
3501  subroutine del2_uv(du, dv, cd, kmd, ntimes, bd, npx, npy, gridstruct, domain)
3502  integer, intent(in):: kmd
3503  integer, intent(in):: ntimes
3504  real, intent(in):: cd
3505  type(fv_grid_bounds_type), intent(IN) :: bd
3506  real, intent(inout), dimension(is:ie,js:je,kmd):: du, dv
3507  integer, intent(IN) :: npx, npy
3508  type(fv_grid_type), intent(IN), target :: gridstruct
3509  type(domain2d), intent(INOUT) :: domain
3510 ! local:
3511  real(kind=R_GRID), pointer, dimension(:,:,:) :: vlon, vlat
3512  real, dimension(is:ie,js:je,kmd):: v1, v2, v3
3513  integer i,j,k
3514 
3515  vlon => gridstruct%vlon
3516  vlat => gridstruct%vlat
3517 
3518 ! transform to 3D Cartesian:
3519 !$OMP parallel do default(none) shared(kmd,is,ie,js,je,v1,v2,v3,du,vlon,dv,vlat)
3520  do k=1,kmd
3521  do j=js,je
3522  do i=is,ie
3523  v1(i,j,k) = du(i,j,k)*vlon(i,j,1) + dv(i,j,k)*vlat(i,j,1)
3524  v2(i,j,k) = du(i,j,k)*vlon(i,j,2) + dv(i,j,k)*vlat(i,j,2)
3525  v3(i,j,k) = du(i,j,k)*vlon(i,j,3) + dv(i,j,k)*vlat(i,j,3)
3526  enddo
3527  enddo
3528  enddo
3529 
3530 ! Filter individual components as scalar:
3531  call del2_scalar( v1(is,js,1), cd, kmd, ntimes, bd, npx, npy, gridstruct, domain )
3532  call del2_scalar( v2(is,js,1), cd, kmd, ntimes, bd, npx, npy, gridstruct, domain )
3533  call del2_scalar( v3(is,js,1), cd, kmd, ntimes, bd, npx, npy, gridstruct, domain )
3534 
3535 ! Convert back to lat-lon components:
3536 !$OMP parallel do default(none) shared(kmd,is,ie,js,je,du,dv,v1,v2,v3,vlon,vlat)
3537  do k=1,kmd
3538  do j=js,je
3539  do i=is,ie
3540  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)
3541  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)
3542  enddo
3543  enddo
3544  enddo
3545 
3546  end subroutine del2_uv
3547 
3549  subroutine del2_scalar(qdt, cd, kmd, nmax, bd, npx, npy, gridstruct, domain)
3550  integer, intent(in):: kmd
3551  integer, intent(in):: nmax
3552  real, intent(in):: cd
3553  type(fv_grid_bounds_type), intent(IN) :: bd
3554  real, intent(inout):: qdt(is:ie,js:je,kmd)
3555  integer, intent(IN) :: npx, npy
3556  type(fv_grid_type), intent(IN), target :: gridstruct
3557  type(domain2d), intent(INOUT) :: domain
3558 ! local:
3559  real:: q(isd:ied,jsd:jed,kmd)
3560  real:: fx(isd:ied+1,jsd:jed), fy(isd:ied,jsd:jed+1)
3561  integer i,j,k, n, nt, ntimes
3562  real :: damp
3563 
3564  real, pointer, dimension(:,:) :: rarea, area
3565  real, pointer, dimension(:,:) :: sina_u, sina_v
3566  real, pointer, dimension(:,:,:) :: sin_sg
3567 
3568  real, pointer, dimension(:,:) :: dx, dy, rdxc, rdyc
3569 
3570  real(kind=R_GRID), pointer :: da_min
3571 
3572  logical, pointer :: nested, sw_corner, se_corner, nw_corner, ne_corner
3573 
3574  area => gridstruct%area
3575  rarea => gridstruct%rarea
3576 
3577  sina_u => gridstruct%sina_u
3578  sina_v => gridstruct%sina_v
3579  sin_sg => gridstruct%sin_sg
3580 
3581  dx => gridstruct%dx
3582  dy => gridstruct%dy
3583  rdxc => gridstruct%rdxc
3584  rdyc => gridstruct%rdyc
3585 
3586  da_min => gridstruct%da_min
3587 
3588  nested => gridstruct%nested
3589  sw_corner => gridstruct%sw_corner
3590  se_corner => gridstruct%se_corner
3591  nw_corner => gridstruct%nw_corner
3592  ne_corner => gridstruct%ne_corner
3593 
3594  ntimes = min(3, nmax)
3595 
3596  damp = cd * da_min
3597 
3598 !$OMP parallel do default(none) shared(is,ie,js,je,kmd,q,qdt)
3599  do k=1,kmd
3600  do j=js,je
3601  do i=is,ie
3602  q(i,j,k) = qdt(i,j,k)
3603  enddo
3604  enddo
3605  enddo
3606  call timing_on('COMM_TOTAL')
3607  call mpp_update_domains(q, domain, complete=.true.)
3608  call timing_off('COMM_TOTAL')
3609 
3610  do n=1,ntimes
3611 
3612  nt = ntimes - n
3613 
3614 !$OMP parallel do default(none) shared(is,ie,js,je,kmd,nt,dy,q,isd,jsd,npx,npy,nested, &
3615 !$OMP bd,sw_corner,se_corner,nw_corner,ne_corner, &
3616 !$OMP sina_u,rdxc,sin_sg,dx,rdyc,sina_v,qdt,damp,rarea) &
3617 !$OMP private(fx, fy)
3618  do k=1,kmd
3619 
3620  if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 1, nested, bd, &
3621  sw_corner, se_corner, nw_corner, ne_corner)
3622  do j=js-nt,je+nt
3623  do i=is-nt,ie+1+nt
3624  fx(i,j) = dy(i,j)*sina_u(i,j)*(q(i-1,j,k)-q(i,j,k))*rdxc(i,j)
3625  enddo
3626  if (is == 1) fx(i,j) = dy(is,j)*(q(is-1,j,k)-q(is,j,k))*rdxc(is,j)* &
3627  0.5*(sin_sg(1,j,1) + sin_sg(0,j,3))
3628  if (ie+1 == npx) fx(i,j) = dy(ie+1,j)*(q(ie,j,k)-q(ie+1,j,k))*rdxc(ie+1,j)* &
3629  0.5*(sin_sg(npx,j,1) + sin_sg(npx-1,j,3))
3630  enddo
3631 
3632  if(nt>0) call copy_corners(q(isd,jsd,k), npx, npy, 2, nested, bd, &
3633  sw_corner, se_corner, nw_corner, ne_corner)
3634  do j=js-nt,je+1+nt
3635  if (j == 1 .OR. j == npy) then
3636  do i=is-nt,ie+nt
3637  fy(i,j) = dx(i,j)*(q(i,j-1,k)-q(i,j,k))*rdyc(i,j) &
3638  *0.5*(sin_sg(i,j,2) + sin_sg(i,j-1,4) )
3639  enddo
3640  else
3641  do i=is-nt,ie+nt
3642  fy(i,j) = dx(i,j)*sina_v(i,j)*(q(i,j-1,k)-q(i,j,k))*rdyc(i,j)
3643  enddo
3644  end if
3645  enddo
3646 
3647  if ( nt==0 ) then
3648  do j=js,je
3649  do i=is,ie
3650  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))
3651  enddo
3652  enddo
3653  else
3654  do j=js-nt,je+nt
3655  do i=is-nt,ie+nt
3656  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))
3657  enddo
3658  enddo
3659  endif
3660  enddo
3661 
3662  enddo
3663 
3664  end subroutine del2_scalar
3665 
3666  subroutine rmse_bias(a, rms, bias, area)
3667  real, intent(in):: a(is:ie,js:je)
3668  real, intent(IN) :: area(isd:ied,jsd:jed)
3669  real, intent(out):: rms, bias
3670  integer:: i,j
3671  real:: total_area
3672 
3673  total_area = 4.*pi*radius**2
3674 
3675  rms = 0.
3676  bias = 0.
3677  do j=js,je
3678  do i=is,ie
3679  bias = bias + area(i,j) * a(i,j)
3680  rms = rms + area(i,j) * a(i,j)**2
3681  enddo
3682  enddo
3683  call mp_reduce_sum(bias)
3684  call mp_reduce_sum(rms)
3685 
3686  bias = bias / total_area
3687  rms = sqrt( rms / total_area )
3688 
3689  end subroutine rmse_bias
3690 
3691 
3692  subroutine corr(a, b, co, area)
3693  real, intent(in):: a(is:ie,js:je), b(is:ie,js:je)
3694  real, intent(in):: area(isd:ied,jsd:jed)
3695  real, intent(out):: co
3696  real:: m_a, m_b, std_a, std_b
3697  integer:: i,j
3698  real:: total_area
3699 
3700  total_area = 4.*pi*radius**2
3701 
3702 ! Compute standard deviation:
3703  call std(a, m_a, std_a, area)
3704  call std(b, m_b, std_b, area)
3705 
3706 ! Compute correlation:
3707  co = 0.
3708  do j=js,je
3709  do i=is,ie
3710  co = co + area(i,j) * (a(i,j)-m_a)*(b(i,j)-m_b)
3711  enddo
3712  enddo
3713  call mp_reduce_sum(co)
3714  co = co / (total_area*std_a*std_b )
3715 
3716  end subroutine corr
3717 
3718  subroutine std(a, mean, stdv, area)
3719  real, intent(in):: a(is:ie,js:je)
3720  real, intent(IN) :: area(isd:ied,jsd:jed)
3721  real, intent(out):: mean, stdv
3722  integer:: i,j
3723  real:: total_area
3724 
3725  total_area = 4.*pi*radius**2
3726 
3727  mean = 0.
3728  do j=js,je
3729  do i=is,ie
3730  mean = mean + area(i,j) * a(i,j)
3731  enddo
3732  enddo
3733  call mp_reduce_sum(mean)
3734  mean = mean / total_area
3735 
3736  stdv = 0.
3737  do j=js,je
3738  do i=is,ie
3739  stdv = stdv + area(i,j) * (a(i,j)-mean)**2
3740  enddo
3741  enddo
3742  call mp_reduce_sum(stdv)
3743  stdv = sqrt( stdv / total_area )
3744 
3745  end subroutine std
3746 
3747 
3748 end module fv_nwp_nudge_mod
real(kind=4), dimension(nobs_max, max_storm) mslp_out
outer ring SLP in mb
Definition: fv_nudge.F90:280
logical do_ps_bias
Definition: fv_nudge.F90:205
real, dimension(:,:), allocatable gz0
Definition: fv_nudge.F90:171
logical module_is_initialized
Definition: fv_nudge.F90:148
real, dimension(:,:,:), allocatable ps_dat
Definition: fv_nudge.F90:168
logical, public do_adiabatic_init
Definition: fv_nudge.F90:138
pure real function, public vicpqd(q)
real, dimension(:), allocatable lon
Definition: fv_nudge.F90:146
subroutine, public get_var3_r4(ncid, var3_name, is, ie, js, je, ks, ke, var3, time_slice)
Definition: sim_nc_mod.F90:246
subroutine breed_srf_winds(time, dt, npz, u_obs, v_obs, ak, bk, ps, phis, delp, ua, va, u_dt, v_dt, pt, q, nwat, zvir, gridstruct)
The subroutine &#39;breed_srf_winds&#39; performs vortex breeding by nudging 1-m winds.
Definition: fv_nudge.F90:2928
character(len=128) analysis_file_first
the first NCEP analysis file to be used for nudging, by default, the first file in the "input_fname_l...
Definition: fv_nudge.F90:176
subroutine compute_slp(isc, iec, jsc, jec, tm, ps, gz, slp)
Definition: fv_nudge.F90:1078
real mask_fac
[0,1] 0: no mask; 1: full strength
Definition: fv_nudge.F90:193
subroutine remap_tq(npz, ak, bk, ps, delp, t, q, kmd, ps0, ta, qa, zvir, ptop)
Definition: fv_nudge.F90:1931
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
subroutine timing_off(blk_name)
The subroutine &#39;timing_off&#39; stops a timer.
Definition: fv_timing.F90:180
The type &#39;fv_grid_type&#39; is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
Definition: fv_arrays.F90:120
integer km
Data z-dimension.
Definition: fv_nudge.F90:144
subroutine del2_uv(du, dv, cd, kmd, ntimes, bd, npx, npy, gridstruct, domain)
The subroutine &#39;del2_uv&#39; filters the wind tendency.
Definition: fv_nudge.F90:3502
integer id_ht_err
Definition: fv_nudge.F90:284
subroutine, public normalize_vect(e)
The subroutine &#39;normalize_vect&#39; makes &#39;e&#39; a unit vector.
subroutine get_ncep_analysis(ps, u, v, t, q, zvir, ts, nfile, fname, bd)
Definition: fv_nudge.F90:1486
real(kind=4), dimension(:,:,:), allocatable gz3
work array
Definition: fv_nudge.F90:170
integer, dimension(:,:), allocatable jdc
Definition: fv_nudge.F90:167
subroutine, public mappm(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord, ptop)
The subroutine &#39;mappm&#39; is a general-purpose routine for remapping one set of vertical levels to anoth...
Definition: fv_mapz.F90:3278
integer, public num_gas
Definition: multi_gases.F90:44
subroutine, public intp_great_circle(beta, p1, p2, x_o, y_o)
real tau_ps
1-day
Definition: fv_nudge.F90:218
subroutine tangent_wind(elon, elat, speed, po, pp, ut, vt)
Definition: fv_nudge.F90:3158
The module &#39;multi_gases&#39; peforms multi constitutents computations.
Definition: multi_gases.F90:25
logical nudge_hght
Definition: fv_nudge.F90:210
subroutine std(a, mean, stdv, area)
Definition: fv_nudge.F90:3719
real, dimension(:,:,:), allocatable s2c
Definition: fv_nudge.F90:166
real(kind=4), dimension(nobs_max, max_storm) time_tc
start time of the track
Definition: fv_nudge.F90:282
real(kind=4), dimension(nobs_max, max_storm) mslp_obs
observed SLP in mb
Definition: fv_nudge.F90:279
real, dimension(:), allocatable bk0
Definition: fv_nudge.F90:145
logical time_varying
Definition: fv_nudge.F90:211
real slp_env
storm environment pressure (pa)
Definition: fv_nudge.F90:252
real pre0_env
critical storm environment pressure (pa) for size computation
Definition: fv_nudge.F90:253
logical strong_mask
Definition: fv_nudge.F90:202
real(kind=r_grid), parameter radius
Definition: fv_nudge.F90:132
integer min_nobs
Definition: fv_nudge.F90:274
real p_norelax
from P_norelax upwards, no nudging
Definition: fv_nudge.F90:184
real(kind=4), dimension(:,:,:,:), allocatable q_dat
Definition: fv_nudge.F90:169
integer, dimension(max_storm) nobs_tc
Definition: fv_nudge.F90:273
integer, parameter, public r_grid
Definition: fv_arrays.F90:35
integer kord_data
Definition: fv_nudge.F90:191
integer jm
Data y-dimension.
Definition: fv_nudge.F90:143
real(kind=4), dimension(nobs_max, max_storm) wind_obs
observed 10-m wind speed (m/s)
Definition: fv_nudge.F90:278
real, parameter del_r
Definition: fv_nudge.F90:262
integer, dimension(:,:), allocatable id2
Definition: fv_nudge.F90:167
integer year_track_data
Definition: fv_nudge.F90:268
pure real function, public virqd(q)
real tau_q
1-day
Definition: fv_nudge.F90:219
The module fv_nwp_nudge contains routines for nudging to input analyses. note This module is currentl...
Definition: fv_nudge.F90:32
subroutine, public get_var1_real(ncid, var1_name, im, var1, var_exist)
Definition: sim_nc_mod.F90:116
integer im
Data x-dimension.
Definition: fv_nudge.F90:142
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
logical nudge_winds
Definition: fv_nudge.F90:208
pure real function, public virq(q)
logical, public t_is_tv
Definition: fv_nudge.F90:195
subroutine get_int_hght(npz, ak, bk, ps, delp, ps0, tv)
Definition: fv_nudge.F90:1895
The module &#39;tp_core&#39; is a collection of routines to support FV transport.
Definition: tp_core.F90:24
subroutine, public fv_nwp_nudge(Time, dt, npx, npy, npz, ps_dt, u_dt, v_dt, t_dt, q_dt, zvir, ptop, ak, bk, ts, ps, delp, ua, va, pt, nwat, q, phis, gridstruct, bd, domain)
Ths subroutine &#39;fv_nwp_nudge&#39; computes and returns time tendencies for nudging to analysis...
Definition: fv_nudge.F90:306
real(kind=4), dimension(:,:,:,:), allocatable v_dat
Definition: fv_nudge.F90:169
subroutine, public get_var1_double(ncid, var1_name, im, var1, var_exist)
The &#39;get_var&#39; subroutines read in variables from netcdf files.
Definition: sim_nc_mod.F90:93
character(len=128) track_file_name
Definition: fv_nudge.F90:188
subroutine remap_coef(agrid)
Definition: fv_nudge.F90:1733
logical use_high_top
Definition: fv_nudge.F90:197
The module &#39;fv_mapz&#39; contains the vertical mapping routines .
Definition: fv_mapz.F90:27
subroutine del2_scalar(qdt, cd, kmd, nmax, bd, npx, npy, gridstruct, domain)
The subroutine &#39;del2_scalar&#39; filters the physics tendency.
Definition: fv_nudge.F90:3550
subroutine, public copy_corners(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
Definition: tp_core.F90:248
logical breed_srf_w
Definition: fv_nudge.F90:245
subroutine, public breed_slp_inline(nstep, dt, npz, ak, bk, phis, pe, pk, peln, pkz, delp, u, v, pt, q, nwat, zvir, gridstruct, ks, domain_local, bd, hydrostatic)
The subroutine &#39;breed_slp_inline&#39; performs vortex breeding by nudging sea level pressure toward singl...
Definition: fv_nudge.F90:2157
subroutine ps_bias_correction(ps_dt, is, ie, js, je, isd, ied, jsd, jed, area)
Definition: fv_nudge.F90:977
subroutine, public fv_nwp_nudge_end
The subroutine &#39;fv_nwp_nudge_end&#39; terminates the nudging module.
Definition: fv_nudge.F90:2085
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
real function, public great_circle_dist(q1, q2, radius)
subroutine pmaxmin(qname, a, imax, jmax, fac)
Definition: fv_nudge.F90:3466
subroutine get_slp_obs(time, nobs, lon_obs, lat_obs, w10, mslp, slp_out, r_out, time_obs, x_o, y_o, w10_o, slp_o, r_vor, p_vor, stime, fact)
Definition: fv_nudge.F90:3187
The module &#39;sim_nc&#39; is a netcdf file reader.
Definition: sim_nc_mod.F90:28
logical print_end_breed
Definition: fv_nudge.F90:212
real function g0_sum(p, ifirst, ilast, jfirst, jlast, mode, reproduce, isd, ied, jsd, jed, area)
Fast version of global sum that is reproduced if the result is rounded to 4-byte. ...
Definition: fv_nudge.F90:1045
real dps_min
maximum PS increment (pa; each step) due to inline breeding
Definition: fv_nudge.F90:163
logical use_pt_inc
Definition: fv_nudge.F90:196
real, dimension(:), allocatable lat
Definition: fv_nudge.F90:146
real, dimension(:), allocatable ak0
Definition: fv_nudge.F90:145
integer, parameter nobs_max
Max # of observations per storm.
Definition: fv_nudge.F90:270
logical nudge_virt
Definition: fv_nudge.F90:209
subroutine remap_uv(npz, ak, bk, ps, delp, u, v, kmd, ps0, u0, v0, ptop)
Definition: fv_nudge.F90:2010
subroutine get_obs(Time, dt, zvir, ak, bk, ps, ts, ps_obs, delp, pt, nwat, q, u_obs, v_obs, t_obs, q_obs, phis, ua, va, u_dt, v_dt, npx, npy, npz, factor, mask, ptop, bd, gridstruct, domain)
Definition: fv_nudge.F90:1095
type(time_type), public fv_time
subroutine, public open_ncfile(iflnm, ncid)
Definition: sim_nc_mod.F90:55
subroutine ncep2fms(sst)
Definition: fv_nudge.F90:1809
subroutine, public vect_cross(e, p1, p2)
The subroutine &#39;vect_cross performs cross products of 3D vectors: e = P1 X P2.
subroutine timing_on(blk_name)
The subroutine &#39;timing_on&#39; starts a timer.
Definition: fv_timing.F90:116
logical nudge_debug
Definition: fv_nudge.F90:204
real p_wvp
cutoff level for specific humidity nudging
Definition: fv_nudge.F90:190
real(kind=4), dimension(nobs_max, max_storm) y_obs
latitude in degrees
Definition: fv_nudge.F90:277
integer, parameter nfile_max
maximum: 20-year analysis data, 4*366*20=29280
Definition: fv_nudge.F90:156
real p_relax
from P_relax upwards, nudging is reduced linearly proportional to pfull/P_relax
Definition: fv_nudge.F90:181
logical print_end_nudge
Definition: fv_nudge.F90:213
subroutine corr(a, b, co, area)
Definition: fv_nudge.F90:3693
real, parameter tm_max
Definition: fv_nudge.F90:254
character(len=128) analysis_file_last
the last NCEP analysis file to be used for nudging by default, the last file in the "input_fname_list...
Definition: fv_nudge.F90:178
integer time_interval
dataset time interval (seconds)
Definition: fv_nudge.F90:153
@ The module &#39;fv_diagnostics&#39; contains routines to compute diagnosic fields.
integer nfile_total
=5 for 1-day (if datasets are 6-hr apart)
Definition: fv_nudge.F90:189
The module &#39;fv_grid_utils&#39; contains routines for setting up and computing grid-related quantities...
logical conserve_mom
Definition: fv_nudge.F90:199
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
logical conserve_hgt
Definition: fv_nudge.F90:200
real nudged_time
seconds
Definition: fv_nudge.F90:264
real function calday(year, month, day, hour, minute, sec)
The function &#39;calday&#39; performs time interpolation:
Definition: fv_nudge.F90:3409
real(kind=4), dimension(:,:,:,:), allocatable t_dat
Definition: fv_nudge.F90:169
subroutine, public fv_nwp_nudge_init(time, axes, npz, zvir, ak, bk, ts, phis, gridstruct, ks, npx, neststruct, bd)
The subroutine &#39;fv_nwp_nudge_init&#39; initializes the nudging module.
Definition: fv_nudge.F90:1257
integer, parameter max_storm
max # of storms to process
Definition: fv_nudge.F90:269
logical nudge_ps
Definition: fv_nudge.F90:206
real(kind=4), dimension(nobs_max, max_storm) rad_out
outer ring radius in meters
Definition: fv_nudge.F90:281
subroutine get_tc_mask(time, mask, agrid)
Definition: fv_nudge.F90:2111
real(kind=4), dimension(nobs_max, max_storm) x_obs
longitude in degrees
Definition: fv_nudge.F90:276
subroutine, public get_ncdim1(ncid, var1_name, im)
Definition: sim_nc_mod.F90:78
subroutine slp_obs_init
Definition: fv_nudge.F90:3279
subroutine breed_srf_w10(time, dt, npz, ak, bk, ps, phis, slp, delp, u, v, gridstruct)
The subroutine &#39;breed_srf_w10&#39; performs vortex breeding by nudging 10-m winds.
Definition: fv_nudge.F90:2649
logical add_bg_wind
Definition: fv_nudge.F90:198
integer kbot_winds
Definition: fv_nudge.F90:236
real tau_winds
6-hr
Definition: fv_nudge.F90:220
character(len=128), dimension(nfile_max) file_names
Definition: fv_nudge.F90:187
integer, dimension(:,:), allocatable id1
Definition: fv_nudge.F90:167
character(len=128) input_fname_list
a file lists the input NCEP analysis data
Definition: fv_nudge.F90:175
logical analysis_time
Definition: fv_nudge.F90:239
logical function leap_year(ny)
Definition: fv_nudge.F90:3443
subroutine rmse_bias(a, rms, bias, area)
Definition: fv_nudge.F90:3667
real(kind=4), dimension(:,:,:,:), allocatable u_dat
Definition: fv_nudge.F90:169
subroutine ps_nudging(dt, factor, npz, ak, bk, ps_obs, mask, tm, ps, phis, delp, ua, va, pt, nwat, q, bd, npx, npy, gridstruct, domain)
Definition: fv_nudge.F90:796
subroutine, public latlon2xyz(p, e, id)
The subroutine &#39;latlon2xyz&#39; maps (lon, lat) to (x,y,z)
subroutine, public close_ncfile(ncid)
Definition: sim_nc_mod.F90:67