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