CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
sfc_nst.f90
1
3
5module sfc_nst
6
7 use machine , only : kind_phys, kp => kind_phys
8 use funcphys , only : fpvs
9 use module_nst_parameters , only : one, zero, half
10 use module_nst_parameters , only : t0k, cp_w, omg_m, omg_sh, sigma_r, solar_time_6am, sst_max
11 use module_nst_parameters , only : ri_c, z_w_max, delz, wd_max, rad2deg, const_rot, tau_min, tw_max
15 !
16 implicit none
17contains
18
27 subroutine sfc_nst_run &
28 ( im, hvap, cp, hfus, jcal, eps, epsm1, rvrdm1, rd, rhw0, & ! --- inputs:
29 pi, tgice, sbc, ps, u1, v1, usfco, vsfco, icplocn2atm, t1, &
30 q1, tref, cm, ch, lseaspray, fm, fm10, &
31 prsl1, prslki, prsik1, prslk1, wet, use_lake_model, xlon, &
32 sinlat, stress, &
33 sfcemis, dlwflx, sfcnsw, rain, timestep, kdt, solhr,xcosz, &
34 wind, flag_iter, flag_guess, nstf_name1, nstf_name4, &
35 nstf_name5, lprnt, ipr, thsfc_loc, &
36 tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, & ! --- input/output:
37 z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain, &
38 qsurf, gflux, cmm, chh, evap, hflx, ep, errmsg, errflg & ! --- outputs:
39 )
40 !
41 ! ===================================================================== !
42 ! description: !
43 ! !
44 ! !
45 ! usage: !
46 ! !
47 ! call sfc_nst !
48 ! inputs: !
49 ! ( im, ps, u1, v1, t1, q1, tref, cm, ch, !
50 ! lseaspray, fm, fm10, !
51 ! prsl1, prslki, wet, use_lake_model, xlon, sinlat, stress, !
52 ! sfcemis, dlwflx, sfcnsw, rain, timestep, kdt,solhr,xcosz, !
53 ! wind, flag_iter, flag_guess, nstf_name1, nstf_name4, !
54 ! nstf_name5, lprnt, ipr, thsfc_loc, !
55 ! input/outputs: !
56 ! tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, !
57 ! z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain, !
58 ! -- outputs:
59 ! qsurf, gflux, cmm, chh, evap, hflx, ep !
60 ! )
61 ! !
62 ! !
63 ! subprogram/functions called: fpvs, density, rhocoef, cool_skin !
64 ! !
65 ! program history log: !
66 ! 2007 -- xu li createad original code !
67 ! 2008 -- s. moorthi adapted to the parallel version !
68 ! may 2009 -- y.-t. hou modified to include input lw surface !
69 ! emissivity from radiation. also replaced the !
70 ! often comfusing combined sw and lw suface !
71 ! flux with separate sfc net sw flux (defined !
72 ! as dn-up) and lw flux. added a program doc block. !
73 ! sep 2009 -- s. moorthi removed rcl and additional reformatting !
74 ! and optimization + made pa as input pressure unit.!
75 ! 2009 -- xu li recreatead the code !
76 ! feb 2010 -- s. moorthi added some changes made to the previous !
77 ! version !
78 ! Jul 2016 -- X. Li, modify the diurnal warming event reset !
79 ! !
80 ! !
81 ! ==================== definition of variables ==================== !
82 ! !
83 ! inputs: size !
84 ! im - integer, horiz dimension 1 !
85 ! ps - real, surface pressure (pa) im !
86 ! u1, v1 - real, u/v component of surface layer wind (m/s) im !
87 ! usfco, vsfco - real, u/v component of surface current (m/s) im !
88 ! icplocn2atm - integer, option to include ocean surface 1 !
89 ! current in the computation of flux !
90 ! t1 - real, surface layer mean temperature ( k ) im !
91 ! q1 - real, surface layer mean specific humidity im !
92 ! tref - real, reference/foundation temperature ( k ) im !
93 ! cm - real, surface exchange coeff for momentum (m/s) im !
94 ! ch - real, surface exchange coeff heat & moisture(m/s) im !
95 ! lseaspray- logical, .t. for parameterization for sea spray 1 !
96 ! fm - real, a stability profile function for momentum im !
97 ! fm10 - real, a stability profile function for momentum im !
98 ! at 10m !
99 ! prsl1 - real, surface layer mean pressure (pa) im !
100 ! prslki - real, im !
101 ! prsik1 - real, im !
102 ! prslk1 - real, im !
103 ! wet - logical, =T if any ocn/lake water (F otherwise) im !
104 ! use_lake_model- logical, =T if flake model is used for lake im !
105 ! icy - logical, =T if any ice im !
106 ! xlon - real, longitude (radians) im !
107 ! sinlat - real, sin of latitude im !
108 ! stress - real, wind stress (n/m**2) im !
109 ! sfcemis - real, sfc lw emissivity (fraction) im !
110 ! dlwflx - real, total sky sfc downward lw flux (w/m**2) im !
111 ! sfcnsw - real, total sky sfc netsw flx into ocean (w/m**2) im !
112 ! rain - real, rainfall rate (kg/m**2/s) im !
113 ! timestep - real, timestep interval (second) 1 !
114 ! kdt - integer, time step counter 1 !
115 ! solhr - real, fcst hour at the end of prev time step 1 !
116 ! xcosz - real, consine of solar zenith angle 1 !
117 ! wind - real, wind speed (m/s) im !
118 ! flag_iter- logical, execution or not im !
119 ! when iter = 1, flag_iter = .true. for all grids im !
120 ! when iter = 2, flag_iter = .true. when wind < 2 im !
121 ! for both land and ocean (when nstf_name1 > 0) im !
122 ! flag_guess-logical, .true.= guess step to get CD et al im !
123 ! when iter = 1, flag_guess = .true. when wind < 2 im !
124 ! when iter = 2, flag_guess = .false. for all grids im !
125 ! nstf_name - integers , NSST related flag parameters 1 !
126 ! nstf_name1 : 0 = NSSTM off 1 !
127 ! 1 = NSSTM on but uncoupled 1 !
128 ! 2 = NSSTM on and coupled 1 !
129 ! nstf_name4 : zsea1 in mm 1 !
130 ! nstf_name5 : zsea2 in mm 1 !
131 ! lprnt - logical, control flag for check print out 1 !
132 ! ipr - integer, grid index for check print out 1 !
133 ! thsfc_loc- logical, flag for reference pressure in theta 1 !
134 ! !
135 ! input/outputs:
136 ! li added for oceanic components
137 ! tskin - real, ocean surface skin temperature ( k ) im !
138 ! tsurf - real, the same as tskin ( k ) but for guess run im !
139 ! xt - real, heat content in dtl im !
140 ! xs - real, salinity content in dtl im !
141 ! xu - real, u-current content in dtl im !
142 ! xv - real, v-current content in dtl im !
143 ! xz - real, dtl thickness im !
144 ! zm - real, mxl thickness im !
145 ! xtts - real, d(xt)/d(ts) im !
146 ! xzts - real, d(xz)/d(ts) im !
147 ! dt_cool - real, sub-layer cooling amount im !
148 ! d_conv - real, thickness of free convection layer (fcl) im !
149 ! z_c - sub-layer cooling thickness im !
150 ! c_0 - coefficient1 to calculate d(tz)/d(ts) im !
151 ! c_d - coefficient2 to calculate d(tz)/d(ts) im !
152 ! w_0 - coefficient3 to calculate d(tz)/d(ts) im !
153 ! w_d - coefficient4 to calculate d(tz)/d(ts) im !
154 ! ifd - real, index to start dtlm run or not im !
155 ! qrain - real, sensible heat flux due to rainfall (watts) im !
156
157 ! outputs: !
158
159 ! qsurf - real, surface air saturation specific humidity im !
160 ! gflux - real, soil heat flux (w/m**2) im !
161 ! cmm - real, im !
162 ! chh - real, im !
163 ! evap - real, evaperation from latent heat flux im !
164 ! hflx - real, sensible heat flux im !
165 ! ep - real, potential evaporation im !
166 ! !
167 ! ===================================================================== !
168
169
170
171 ! --- inputs:
172 integer, intent(in) :: im, kdt, ipr, nstf_name1, nstf_name4, nstf_name5
173 integer, intent(in) :: icplocn2atm
174
175 real (kind=kind_phys), intent(in) :: hvap, cp, hfus, jcal, eps, &
176 epsm1, rvrdm1, rd, rhw0, sbc, pi, tgice
177 real (kind=kind_phys), dimension(:), intent(in) :: ps, u1, v1, &
178 usfco, vsfco, t1, q1, cm, ch, fm, fm10, &
179 prsl1, prslki, prsik1, prslk1, xlon, xcosz, &
180 sinlat, stress, sfcemis, dlwflx, sfcnsw, rain, wind
181 real (kind=kind_phys), dimension(:), intent(in), optional :: &
182 tref
183 real (kind=kind_phys), intent(in) :: timestep
184 real (kind=kind_phys), intent(in) :: solhr
185
186 ! For sea spray effect
187 logical, intent(in) :: lseaspray
188 !
189 logical, dimension(:), intent(in) :: flag_iter, flag_guess, wet
190 integer, dimension(:), intent(in) :: use_lake_model
191 logical, intent(in) :: lprnt
192 logical, intent(in) :: thsfc_loc
193
194 ! --- input/outputs:
195 ! control variables of dtl system (5+2) and sl (2) and coefficients for d(tz)/d(ts) calculation
196 real (kind=kind_phys), dimension(:), intent(inout) :: tskin, &
197 tsurf
198 real (kind=kind_phys), dimension(:), intent(inout), optional :: &
199 xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, &
200 z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain
201
202 ! --- outputs:
203 real (kind=kind_phys), dimension(:), intent(inout) :: qsurf, gflux, cmm, chh, evap, hflx, ep
204
205 character(len=*), intent(out) :: errmsg
206 integer, intent(out) :: errflg
207
208 !
209 ! locals
210 !
211 integer :: k,i
212 !
213 real (kind=kind_phys), dimension(im) :: q0, qss, rch, rho_a, theta1, tv1, wndmag
214
215 real(kind=kind_phys) :: elocp,tem,cpinv,hvapi
216 !
217 ! nstm related prognostic fields
218 !
219 logical :: flag(im)
220 real (kind=kind_phys), dimension(im) :: xt_old, xs_old, xu_old, xv_old, xz_old, &
221 zm_old,xtts_old, xzts_old, ifd_old, tref_old, tskin_old, dt_cool_old,z_c_old
222
223 real(kind=kind_phys) :: ulwflx(im), nswsfc(im)
224 ! real(kind=kind_phys) rig(im),
225 ! & ulwflx(im),dlwflx(im),
226 ! & slrad(im),nswsfc(im)
227 real(kind=kind_phys) :: alpha,beta,rho_w,f_nsol,sss,sep, cosa,sina,taux,tauy, &
228 grav,dz,t0,ttop0,ttop
229
230 real(kind=kind_phys) :: le,fc,dwat,dtmp,wetc,alfac,ustar_a,rich
231 real(kind=kind_phys) :: rnl_ts,hs_ts,hl_ts,rf_ts,q_ts
232 real(kind=kind_phys) :: fw,q_warm
233 real(kind=kind_phys) :: t12,alon,tsea,sstc,dta,dtz
234 real(kind=kind_phys) :: zsea1,zsea2,soltim
235 logical :: do_nst
236 !
237 ! parameters for sea spray effect
238 !
239 real (kind=kind_phys) :: f10m, u10m, v10m, ws10, ru10, qss1, &
240 bb1, hflxs, evaps, ptem
241 !
242 ! real (kind=kind_phys), parameter :: alps=0.5, bets=0.5, gams=0.1,
243 ! real (kind=kind_phys), parameter :: alps=0.5, bets=0.5, gams=0.0,
244 ! real (kind=kind_phys), parameter :: alps=1.0, bets=1.0, gams=0.2,
245 real (kind=kind_phys), parameter :: alps=0.75,bets=0.75,gams=0.15, &
246 ws10cr=30., conlf=7.2e-9, consf=6.4e-8
247 real (kind=kind_phys) :: windrel
248 !
249 !======================================================================================================
250 ! Initialize CCPP error handling variables
251 errmsg = ''
252 errflg = 0
253
254 if (nstf_name1 == 0) return ! No NSST model used
255
256 cpinv = one/cp
257 hvapi = one/hvap
258 elocp = hvap/cp
259
260 sss = 34.0_kp ! temporarily, when sea surface salinity data is not ready
261 !
262 ! flag for open water and where the iteration is on
263 !
264 do_nst = .false.
265 do i = 1, im
266 ! flag(i) = wet(i) .and. .not.icy(i) .and. flag_iter(i)
267 flag(i) = wet(i) .and. flag_iter(i) .and. use_lake_model(i)/=1
268 do_nst = do_nst .or. flag(i)
269 enddo
270 if (.not. do_nst) return
271 !
272 ! save nst-related prognostic fields for guess run
273 !
274 do i=1, im
275 ! if(wet(i) .and. .not.icy(i) .and. flag_guess(i)) then
276 if(wet(i) .and. flag_guess(i) .and. use_lake_model(i)/=1) then
277 xt_old(i) = xt(i)
278 xs_old(i) = xs(i)
279 xu_old(i) = xu(i)
280 xv_old(i) = xv(i)
281 xz_old(i) = xz(i)
282 zm_old(i) = zm(i)
283 xtts_old(i) = xtts(i)
284 xzts_old(i) = xzts(i)
285 ifd_old(i) = ifd(i)
286 tskin_old(i) = tskin(i)
287 dt_cool_old(i) = dt_cool(i)
288 z_c_old(i) = z_c(i)
289 endif
290 enddo
291
292
293 ! --- ... initialize variables. all units are m.k.s. unless specified.
294 ! ps is in pascals, wind is wind speed, theta1 is surface air
295 ! estimated from level 1 temperature, rho_a is air density and
296 ! qss is saturation specific humidity at the water surface
297 !!
298 do i = 1, im
299 if ( flag(i) ) then
300
301 nswsfc(i) = sfcnsw(i) ! net solar radiation at the air-sea surface (positive downward)
302 wndmag(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i))
303
304 q0(i) = max(q1(i), 1.0e-8_kp)
305
306 if(thsfc_loc) then ! Use local potential temperature
307 theta1(i) = t1(i) * prslki(i)
308 else ! Use potential temperature referenced to 1000 hPa
309 theta1(i) = t1(i) / prslk1(i) ! potential temperature at the middle of lowest model layer
310 endif
311
312 tv1(i) = t1(i) * (one + rvrdm1*q0(i))
313 rho_a(i) = prsl1(i) / (rd*tv1(i))
314 qss(i) = fpvs(tsurf(i)) ! pa
315 qss(i) = eps*qss(i) / (ps(i) + epsm1*qss(i)) ! pa
316 !
317 evap(i) = zero
318 hflx(i) = zero
319 gflux(i) = zero
320 ep(i) = zero
321
322 ! --- ... rcp = rho cp ch v
323
324 if (icplocn2atm ==0) then
325 rch(i) = rho_a(i) * cp * ch(i) * wind(i)
326 cmm(i) = cm(i) * wind(i)
327 chh(i) = rho_a(i) * ch(i) * wind(i)
328 else if (icplocn2atm ==1) then
329 windrel= sqrt( (u1(i)-usfco(i))**2 + (v1(i)-vsfco(i))**2 )
330 rch(i) = rho_a(i) * cp * ch(i) * windrel
331 cmm(i) = cm(i) * windrel
332 chh(i) = rho_a(i) * ch(i) * windrel
333 endif
334
336 ! at previous time step
337 evap(i) = elocp * rch(i) * (qss(i) - q0(i))
338 qsurf(i) = qss(i)
339
340 if(thsfc_loc) then ! Use local potential temperature
341 hflx(i) = rch(i) * (tsurf(i) - theta1(i))
342 else ! Use potential temperature referenced to 1000 hPa
343 hflx(i) = rch(i) * (tsurf(i)/prsik1(i) - theta1(i))
344 endif
345
346 ! if (lprnt .and. i == ipr) print *,' tskin=',tskin(i),' theta1=',
347 ! & theta1(i),' hflx=',hflx(i),' t1=',t1(i),'prslki=',prslki(i)
348 ! &,' tsurf=',tsurf(i)
349 endif
350 enddo
351
352 ! run nst model: dtm + slm
353 !
354 zsea1 = 0.001_kp*real(nstf_name4)
355 zsea2 = 0.001_kp*real(nstf_name5)
356
360 do i = 1, im
361 if ( flag(i) ) then
362 tsea = tsurf(i)
363 t12 = tsea*tsea
364 ulwflx(i) = sfcemis(i) * sbc * t12 * t12
365 alon = xlon(i)*rad2deg
366 grav = grv(sinlat(i))
367 soltim = mod(alon/15.0_kp + solhr, 24.0_kp)*3600.0_kp
368 call density(tsea,sss,rho_w) ! sea water density
369 call rhocoef(tsea,sss,rho_w,alpha,beta) ! alpha & beta
370 !
372 !
373 le = (2.501_kp-0.00237_kp*tsea)*1.0e6_kp
374 dwat = 2.11e-5_kp*(t1(i)/t0k)**1.94_kp ! water vapor diffusivity
375 dtmp = (one+3.309e-3_kp*(t1(i)-t0k)-1.44e-6_kp*(t1(i)-t0k) &
376 * (t1(i)-t0k))*0.02411_kp/(rho_a(i)*cp) ! heat diffusivity
377 wetc = 622.0_kp*le*qss(i)/(rd*t1(i)*t1(i))
378 alfac = one / (one + (wetc*le*dwat)/(cp*dtmp)) ! wet bulb factor
379 tem = (1.0e3_kp * rain(i) / rho_w) * alfac * cp_w
380 qrain(i) = tem * (tsea-t1(i)+1.0e3_kp*(qss(i)-q0(i))*le/cp)
381
383
384 f_nsol = hflx(i) + evap(i) + ulwflx(i) - dlwflx(i) + omg_sh*qrain(i)
385
386 ! if (lprnt .and. i == ipr) print *,' f_nsol=',f_nsol,' hflx=',
387 ! &hflx(i),' evap=',evap(i),' ulwflx=',ulwflx(i),' dlwflx=',dlwflx(i)
388 ! &,' omg_sh=',omg_sh,' qrain=',qrain(i)
389
390 sep = sss*(evap(i)/le-rain(i))/rho_w
391 ustar_a = sqrt(stress(i)/rho_a(i)) ! air friction velocity
392 !
393 ! sensitivities of heat flux components to ts
394 !
395 rnl_ts = 4.0_kp*sfcemis(i)*sbc*tsea*tsea*tsea ! d(rnl)/d(ts)
396 hs_ts = rch(i)
397 hl_ts = rch(i)*elocp*eps*hvap*qss(i)/(rd*t12)
398 rf_ts = tem * (one+rch(i)*hl_ts)
399 q_ts = rnl_ts + hs_ts + hl_ts + omg_sh*rf_ts
400 !
403 ! & calculate c_0, c_d
404 !
405 call cool_skin(ustar_a,f_nsol,nswsfc(i),evap(i),sss,alpha,beta, &
406 rho_w,rho_a(i),tsea,q_ts,hl_ts,grav,le, &
407 dt_cool(i),z_c(i),c_0(i),c_d(i))
408
409 tem = one / wndmag(i)
410 cosa = u1(i)*tem
411 sina = v1(i)*tem
412 taux = max(stress(i),tau_min)*cosa
413 tauy = max(stress(i),tau_min)*sina
414 fc = const_rot*sinlat(i)
415 !
416 ! Run DTM-1p system.
417 !
418 if ( (soltim > solar_time_6am .and. ifd(i) == zero) ) then
419 else
420 ifd(i) = one
421 !
422 ! calculate fcl thickness with current forcing and previous time's profile
423 !
424 ! if (lprnt .and. i == ipr) print *,' beg xz=',xz(i)
425
427 if ( f_nsol > zero .and. xt(i) > zero ) then
428 call convdepth(kdt,timestep,nswsfc(i),f_nsol,sss,sep,rho_w, &
429 alpha,beta,xt(i),xs(i),xz(i),d_conv(i))
430 else
431 d_conv(i) = zero
432 endif
433
434 ! if (lprnt .and. i == ipr) print *,' beg xz1=',xz(i)
435 !
436 ! determine rich: wind speed dependent (right now)
437 !
438 ! if ( wind(i) < 1.0 ) then
439 ! rich = 0.25 + 0.03*wind(i)
440 ! elseif ( wind(i) >= 1.0 .and. wind(i) < 1.5 ) then
441 ! rich = 0.25 + 0.1*wind(i)
442 ! elseif ( wind(i) >= 1.5 .and. wind(i) < 6.0 ) then
443 ! rich = 0.25 + 0.6*wind(i)
444 ! elseif ( wind(i) >= 6.0 ) then
445 ! rich = 0.25 + min(0.8*wind(i),0.50)
446 ! endif
447
448 rich = ri_c
449
451 call dtm_1p(kdt,timestep,rich,taux,tauy,nswsfc(i), &
452 f_nsol,sss,sep,q_ts,hl_ts,rho_w,alpha,beta,alon, &
453 sinlat(i),soltim,grav,le,d_conv(i), &
454 xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
455
456 ! if (lprnt .and. i == ipr) print *,' beg xz2=',xz(i)
457
458 ! apply mda
459 if ( xt(i) > zero ) then
462 call dtm_1p_mda(xt(i),xtts(i),xz(i),xzts(i))
463 if ( xz(i) >= z_w_max ) then
466 call dtl_reset(xt(i),xs(i),xu(i),xv(i),xz(i),xtts(i), xzts(i))
467
468 ! if (lprnt .and. i == ipr) print *,' beg xz3=',xz(i),' z_w_max='
469 ! &,z_w_max
470 endif
471
472 ! apply fca
473 if ( d_conv(i) > zero ) then
478 call dtm_1p_fca(d_conv(i),xt(i),xtts(i),xz(i),xzts(i))
479 if ( xz(i) >= z_w_max ) then
480 call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
481 endif
482 endif
483
484 ! if (lprnt .and. i == ipr) print *,' beg xz4=',xz(i)
485
486 ! apply tla
487 dz = min(xz(i),max(d_conv(i),delz))
488 !
492 call sw_ps_9b(delz,fw)
493 q_warm = fw*nswsfc(i)-f_nsol !total heat absorbed in warm layer
494
497 if ( q_warm > zero ) then
498 call cal_ttop(kdt,timestep,q_warm,rho_w,dz, xt(i),xz(i),ttop0)
499
500 ! if (lprnt .and. i == ipr) print *,' d_conv=',d_conv(i),' delz=',
501 ! &delz,' kdt=',kdt,' timestep=',timestep,' nswsfc=',nswsfc(i),
502 ! &' f_nsol=',f_nsol,' rho_w=',rho_w,' dz=',dz,' xt=',xt(i),
503 ! &' xz=',xz(i),' qrain=',qrain(i)
504
505 ttop = ((xt(i)+xt(i))/xz(i))*(one-dz/((xz(i)+xz(i))))
506
507 ! if (lprnt .and. i == ipr) print *,' beg xz4a=',xz(i)
508 ! &,' ttop=',ttop,' ttop0=',ttop0,' xt=',xt(i),' dz=',dz
509 ! &,' xznew=',(xt(i)+sqrt(xt(i)*(xt(i)-dz*ttop0)))/ttop0
510
512 if ( ttop > ttop0 ) then
513 call dtm_1p_tla(dz,ttop0,xt(i),xtts(i),xz(i),xzts(i))
514
515 ! if (lprnt .and. i == ipr) print *,' beg xz4b=',xz(i),'z_w_max=',
516 ! &z_w_max
517 if ( xz(i) >= z_w_max ) then
518 call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
519 endif
520 endif
521 endif ! if ( q_warm > 0.0 ) then
522
523 ! if (lprnt .and. i == ipr) print *,' beg xz5=',xz(i)
524
525 ! apply mwa
527 t0 = (xt(i)+xt(i))/xz(i)
528 if ( t0 > tw_max ) then
529 call dtm_1p_mwa(xt(i),xtts(i),xz(i),xzts(i))
530 if ( xz(i) >= z_w_max ) then
531 call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
532 endif
533 endif
534
535 ! if (lprnt .and. i == ipr) print *,' beg xz6=',xz(i)
536
537 ! apply mta
539 sstc = tref(i) + (xt(i)+xt(i))/xz(i) - dt_cool(i)
540
541 if ( sstc > sst_max ) then
542 dta = sstc - sst_max
543 call dtm_1p_mta(dta,xt(i),xtts(i),xz(i),xzts(i))
544 ! write(*,'(a,f3.0,7f8.3)') 'mta, sstc,dta :',islimsk(i),
545 ! & sstc,dta,tref(i),xt(i),xz(i),2.0*xt(i)/xz(i),dt_cool(i)
546 if ( xz(i) >= z_w_max ) then
547 call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
548 endif
549 endif
550 !
551 endif ! if ( xt(i) > 0.0 ) then
552 ! reset dtl at midnight and when solar zenith angle > 89.994 degree
553 if ( abs(soltim) < 2.0_kp*timestep ) then
554 call dtl_reset (xt(i),xs(i),xu(i),xv(i),xz(i),xzts(i),xtts(i))
555 endif
556
557 endif ! if (solar_time > solar_time_6am .and. ifd(i) == 0.0 ) then: too late to start the first day
558
559 ! if (lprnt .and. i == ipr) print *,' beg xz7=',xz(i)
560
561 ! update tsurf (when flag(i) .eqv. .true. )
563 call get_dtzm_point(xt(i),xz(i),dt_cool(i),z_c(i), zsea1,zsea2,dtz)
564 tsurf(i) = max(tgice, tref(i) + dtz )
565
566 ! if (lprnt .and. i == ipr) print *,' tsurf=',tsurf(i),' tref=',
567 ! &tref(i),' xz=',xz(i),' dt_cool=',dt_cool(i)
568
570 if ( xt(i) > zero ) then
571 call cal_w(kdt,xz(i),xt(i),xzts(i),xtts(i),w_0(i),w_d(i))
572 else
573 w_0(i) = zero
574 w_d(i) = zero
575 endif
576
577 ! if ( xt(i) > 0.0 ) then
578 ! rig(i) = grav*xz(i)*xz(i)*(alpha*xt(i)-beta*xs(i))
579 ! & /(2.0*(xu(i)*xu(i)+xv(i)*xv(i)))
580 ! else
581 ! rig(i) = 0.25
582 ! endif
583
584 ! qrain(i) = rig(i)
585 zm(i) = wind(i)
586
587 endif
588 enddo
589
590 ! restore nst-related prognostic fields for guess run
591 do i=1, im
592 ! if (wet(i) .and. .not.icy(i)) then
593 if (wet(i) .and. use_lake_model(i)/=1) then
594 if (flag_guess(i)) then ! when it is guess of
595 xt(i) = xt_old(i)
596 xs(i) = xs_old(i)
597 xu(i) = xu_old(i)
598 xv(i) = xv_old(i)
599 xz(i) = xz_old(i)
600 zm(i) = zm_old(i)
601 xtts(i) = xtts_old(i)
602 xzts(i) = xzts_old(i)
603 ifd(i) = ifd_old(i)
604 tskin(i) = tskin_old(i)
605 dt_cool(i) = dt_cool_old(i)
606 z_c(i) = z_c_old(i)
607 else
608 !
609 ! update tskin when coupled and not guess run
610 ! (all other NSST variables have been updated in this case)
611 !
612 if ( nstf_name1 > 1 ) then
613 tskin(i) = tsurf(i)
614 endif ! if nstf_name1 > 1 then
615 endif ! if flag_guess(i) then
616 endif ! if wet(i) .and. .not.icy(i) then
617 enddo
618
619 ! if (lprnt .and. i == ipr) print *,' beg xz8=',xz(i)
620
621 if ( nstf_name1 > 1 ) then
624 do i = 1, im
625 if ( flag(i) ) then
626 qss(i) = fpvs( tskin(i) )
627 qss(i) = eps*qss(i) / (ps(i) + epsm1*qss(i))
628 qsurf(i) = qss(i)
629 evap(i) = elocp*rch(i) * (qss(i) - q0(i))
630
631 if(thsfc_loc) then ! Use local potential temperature
632 hflx(i) = rch(i) * (tskin(i) - theta1(i))
633 else ! Use potential temperature referenced to 1000 hPa
634 hflx(i) = rch(i) * (tskin(i)/prsik1(i) - theta1(i))
635 endif
636
637 endif
638 enddo
639 endif ! if ( nstf_name1 > 1 ) then
640 !
642 !
643 do i=1,im
644 if(lseaspray .and. flag(i)) then
645 f10m = fm10(i) / fm(i)
646 u10m = f10m * u1(i)
647 v10m = f10m * v1(i)
648 ws10 = sqrt(u10m*u10m + v10m*v10m)
649 ws10 = max(ws10,1.)
650 ws10 = min(ws10,ws10cr)
651 tem = .015 * ws10 * ws10
652 ru10 = 1. - .087 * log(10./tem)
653 qss1 = fpvs(t1(i))
654 qss1 = eps * qss1 / (prsl1(i) + epsm1 * qss1)
655 tem = rd * cp * t1(i) * t1(i)
656 tem = 1. + eps * hvap * hvap * qss1 / tem
657 bb1 = 1. / tem
658 evaps = conlf * (ws10**5.4) * ru10 * bb1
659 evaps = evaps * rho_a(i) * hvap * (qss1 - q0(i))
660 evap(i) = evap(i) + alps * evaps
661 hflxs = consf * (ws10**3.4) * ru10
662 hflxs = hflxs * rho_a(i) * cp * (tskin(i) - t1(i))
663 ptem = alps - gams
664 hflx(i) = hflx(i) + bets * hflxs - ptem * evaps
665 endif
666 enddo
667 !
668 do i=1,im
669 if ( flag(i) ) then
670 tem = one / rho_a(i)
671 hflx(i) = hflx(i) * tem * cpinv
672 evap(i) = evap(i) * tem * hvapi
673 endif
674 enddo
675 !
676 ! if (lprnt) print *,' tskin=',tskin(ipr)
677
678 return
679 end subroutine sfc_nst_run
681end module sfc_nst
subroutine, public convdepth(kdt, timestep, i0, q, sss, sep, rho, alpha, beta, xt, xs, xz, d_conv)
This subroutine calculates depth for convective adjustment.
subroutine, public dtm_1p_mwa(xt, xtts, xz, xzts)
This subroutine applies maximum warming adjustment (mwa).
subroutine, public density(t, s, rho)
This subroutine computes sea water density.
subroutine, public dtm_1p_fca(d_conv, xt, xtts, xz, xzts)
This subroutine applies free convection adjustment(fca).
subroutine, public get_dtzm_point(xt, xz, dt_cool, zc, z1, z2, dtm)
This subroutine computes dtm (the mean of ).
subroutine, public cal_ttop(kdt, timestep, q_warm, rho, dz, xt, xz, ttop)
This subroutine calculates the diurnal warming amount at the top layer with thickness of delz.
subroutine, public dtm_1p_mta(dta, xt, xtts, xz, xzts)
This subroutine applies maximum temperature adjustment (mta).
subroutine, public dtm_1p(kdt, timestep, rich, tox, toy, i0, q, sss, sep, q_ts, hl_ts, rho, alpha, beta, alon, sinlat, soltim, grav, le, d_conv, xt, xs, xu, xv, xz, xzts, xtts)
This subroutine contains the module of diurnal thermocline layer model.
subroutine sfc_nst_run(im, hvap, cp, hfus, jcal, eps, epsm1, rvrdm1, rd, rhw0, pi, tgice, sbc, ps, u1, v1, usfco, vsfco, icplocn2atm, t1, q1, tref, cm, ch, lseaspray, fm, fm10, prsl1, prslki, prsik1, prslk1, wet, use_lake_model, xlon, sinlat, stress, sfcemis, dlwflx, sfcnsw, rain, timestep, kdt, solhr, xcosz, wind, flag_iter, flag_guess, nstf_name1, nstf_name4, nstf_name5, lprnt, ipr, thsfc_loc, tskin, tsurf, xt, xs, xu, xv, xz, zm, xtts, xzts, dt_cool, z_c, c_0, c_d, w_0, w_d, d_conv, ifd, qrain, qsurf, gflux, cmm, chh, evap, hflx, ep, errmsg, errflg)
Definition sfc_nst.f90:40
subroutine, public dtm_1p_mda(xt, xtts, xz, xzts)
This subroutine applies minimum depth adjustment (xz adjustment).
subroutine, public dtl_reset(xt, xs, xu, xv, xz, xzts, xtts)
This subroutine resets the value of xt,xs,xu,xv,xz,xtts,xzts.
subroutine, public rhocoef(t, s, rhoref, alpha, beta)
This subroutine computes thermal expansion coefficient (alpha) and saline contraction coefficient (be...
subroutine, public dtm_1p_tla(dz, te, xt, xtts, xz, xzts)
This subroutine applies top layer adjustment (tla).
subroutine, public cal_w(kdt, xz, xt, xzts, xtts, w_0, w_d)
This subroutine computes coefficients (w_0 and w_d) to calculate d(tw)/d(ts).
subroutine, public cool_skin(ustar_a, f_nsol, f_sol_0, evap, sss, alpha, beta, rho_w, rho_a, ts, q_ts, hl_ts, grav, le, deltat_c, z_c, c_0, c_d)
This subroutine contains the upper ocean cool-skin parameterization (Fairall et al,...
real(kind_phys) function, public grv(x)
This module contains the diurnal thermocline layer model (DTM) of the GFS NSST scheme.
This module contains the CCPP-compliant GFS near-surface sea temperature scheme.
Definition sfc_nst.f90:5