CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_nst_model.f90
1
4
9
13 !
14 ! the module of diurnal thermocline layer model
15 !
16 use machine , only : kind_phys
17 use module_nst_parameters , only : z_w_max, z_w_min, z_w_ini, eps_z_w, eps_conv
18 use module_nst_parameters , only : eps_sfs, niter_z_w, niter_conv, niter_sfs, ri_c
19 use module_nst_parameters , only : ri_g, omg_m, omg_sh, kw => tc_w, visw, t0k, cp_w
20 use module_nst_parameters , only : z_c_max, z_c_ini, ustar_a_min, delz, exp_const
21 use module_nst_parameters , only : rad2deg, const_rot, tw_max, sst_max
22 use module_nst_parameters , only : zero, one
24
25 implicit none
26
27 private
28
31
32contains
33
36 subroutine dtm_1p(kdt,timestep,rich,tox,toy,i0,q,sss,sep,q_ts,hl_ts,rho, &
37 alpha,beta,alon,sinlat,soltim,grav,le,d_conv, &
38 xt,xs,xu,xv,xz,xzts,xtts)
39
40 integer, intent(in) :: kdt
41 real(kind=kind_phys), intent(in) :: timestep,rich,tox,toy,i0,q,sss,sep,q_ts,&
42 hl_ts,rho,alpha,beta,alon,sinlat,soltim,&
43 grav,le,d_conv
44 real(kind=kind_phys), intent(inout) :: xt,xs,xu,xv,xz,xzts,xtts
45 ! local variables
46
47 !
48 ! input variables
49 !
50 ! timestep: integration time step in seconds
51 ! rich : critical ri (flow dependent)
52 ! tox : x wind stress (n*m^-2 or kg/m/s^2)
53 ! toy : y wind stress (n*m^-2 or kg/m/s^2)
54 ! i0 : solar radiation flux at the surface (wm^-2)
55 ! q : non-solar heat flux at the surface (wm^-2)
56 ! sss : salinity (ppt)
57 ! sep : sr(e-p) (ppt*m/s)
58 ! q_ts : d(q)/d(ts) : q = the sum of non-solar heat fluxes
59 ! hl_ts : d(hl)/d(ts)
60 ! rho : sea water density (kg*m^-3)
61 ! alpha : thermal expansion coefficient (1/k)
62 ! beta : saline contraction coefficient (1/ppt)
63 ! sinlat : sine (lat)
64 ! grav : gravity accelleration
65 ! le : le=(2.501-.00237*tsea)*1e6
66 ! d-conv : fcl thickness
67 !
68 ! inout variables
69 !
70 ! xt : dtl heat content (m*k)
71 ! xs : dtl salinity content (m*ppt)
72 ! xu : dtl x current content (m*m/s)
73 ! xv : dtl y current content (m*m/s)
74 ! xz : dtl thickness (m)
75 ! xzts : d(xz)/d(ts) (m/k )
76 ! xtts : d(xt)/d(ts) (m)
77 !
78 ! logical lprnt
79
80 ! if (lprnt) print *,' first xt=',xt
81 if ( xt <= zero ) then ! dtl doesn't exist yet
82 call dtm_onset(kdt,timestep,rich,tox,toy,i0,q,sss,sep,q_ts,hl_ts,rho,alpha, &
83 beta,alon,sinlat,soltim,grav,le,xt,xs,xu,xv,xz,xzts,xtts)
84 elseif ( xt > zero ) then ! dtl already exists
85 !
86 ! forward the system one time step
87 !
88 call eulerm(kdt,timestep,rich,tox,toy,i0,q,sss,sep,q_ts,hl_ts,rho,alpha, &
89 beta,alon,sinlat,soltim,grav,le,d_conv, &
90 xt,xs,xu,xv,xz,xzts,xtts)
91 endif ! if ( xt == 0 ) then
92
93 end subroutine dtm_1p
94
97 subroutine eulerm(kdt,timestep,rich,tox,toy,i0,q,sss,sep,q_ts,hl_ts,rho,alpha, &
98 beta,alon,sinlat,soltim,grav,le,d_conv, &
99 xt,xs,xu,xv,xz,xzts,xtts)
100
101 !
102 ! subroutine eulerm: integrate one time step with modified euler method
103 !
104 integer, intent(in) :: kdt
105 real(kind=kind_phys), intent(in) :: timestep,rich,tox,toy,i0,q,sss,sep,q_ts, &
106 hl_ts,rho,alpha,beta,alon,sinlat,soltim, &
107 grav,le,d_conv
108 real(kind=kind_phys), intent(inout) :: xt,xs,xu,xv,xz,xzts,xtts
109 ! local variables
110 real(kind=kind_phys) :: xt0,xs0,xu0,xv0,xz0,xzts0,xtts0
111 real(kind=kind_phys) :: fw,aw,q_warm
112 real(kind=kind_phys) :: xt1,xs1,xu1,xv1,xz1,xzts1,xtts1
113 real(kind=kind_phys) :: xt2,xs2,xu2,xv2,xz2,xzts2,xtts2
114 real(kind=kind_phys) :: dzw,drho,fc
115 real(kind=kind_phys) :: alat,speed
116 ! logical lprnt
117
118 !
119 ! input variables
120 !
121 ! timestep: integration time step in seconds
122 ! rich : critial ri (flow/mass dependent)
123 ! tox : x wind stress (n*m^-2 or kg/m/s^2)
124 ! toy : y wind stress (n*m^-2 or kg/m/s^2)
125 ! i0 : solar radiation flux at the surface (wm^-2)
126 ! q : non-solar heat flux at the surface (wm^-2)
127 ! sss : salinity (ppt)
128 ! sep : sr(e-p) (ppt*m/s)
129 ! q_ts : d(q)/d(ts) : q = the sum of non-solar heat fluxes
130 ! hl_ts : d(hl)/d(ts)
131 ! rho : sea water density (kg*m^-3)
132 ! alpha : thermal expansion coefficient (1/k)
133 ! beta : saline contraction coefficient (1/ppt)
134 ! alon : longitude (deg)
135 ! sinlat : sine (lat)
136 ! soltim : solar time
137 ! grav : gravity accelleration
138 ! le : le=(2.501-.00237*tsea)*1e6
139 ! d_conv : fcl thickness (m)
140 !
141 ! inout variables
142 !
143 ! xt : dtl heat content (m*k)
144 ! xs : dtl salinity content (m*ppt)
145 ! xu : dtl x current content (m*m/s)
146 ! xv : dtl y current content (m*m/s)
147 ! xz : dtl thickness (m)
148 ! xzts : d(xz)/d(ts) (m/k )
149 ! xtts : d(xt)/d(ts) (m)
150
151 xt0 = xt
152 xs0 = xs
153 xu0 = xu
154 xv0 = xv
155 xz0 = xz
156 xtts0 = xtts
157 xzts0 = xzts
158 speed = max(1.0e-8, xu0*xu0+xv0*xv0)
159
160 alat = asin(sinlat)*rad2deg
161
162 fc = const_rot*sinlat
163
164 call sw_ps_9b(xz0,fw)
165
166 q_warm = fw*i0-q !total heat abs in warm layer
167
168 call sw_ps_9b_aw(xz0,aw)
169
170 drho = -alpha*q_warm/(rho*cp_w) + omg_m*beta*sep
171
172 ! dzw = xz0*(tox*xu0+toy*xv0) / (rho*(xu0*xu0+xv0*xv0)) &
173 ! + xz0*xz0*xz0*drho*grav / (4.0*rich*(xu0*xu0+xv0*xv0))
174 dzw = xz0 * ((tox*xu0+toy*xv0) / (rho*speed) &
175 + xz0*xz0*drho*grav / (4.0*rich*speed))
176
177 xt1 = xt0 + timestep*q_warm/(rho*cp_w)
178 xs1 = xs0 + timestep*sep
179 xu1 = xu0 + timestep*(fc*xv0+tox/rho)
180 xv1 = xv0 + timestep*(-fc*xu0+toy/rho)
181 xz1 = xz0 + timestep*dzw
182
183 ! if (lprnt) print *,' xt1=',xt1,' xz1=',xz1,' xz0=',xz0,' dzw=',dzw, &
184 ! 'timestep=',timestep,tox,toy,xu0,xv0,rho,drho,grav,rich
185
186 if ( xt1 <= zero .or. xz1 <= zero .or. xz1 > z_w_max ) then
187 call dtl_reset(xt,xs,xu,xv,xz,xzts,xtts)
188 return
189 endif
190
191 ! call dtm_1p_zwa(kdt,timestep,i0,q,rho,d_conv,xt1,xs1,xu1,xv1,xz1,tr_mda,tr_fca,tr_tla,tr_mwa)
192
193 xzts1 = xzts0 + timestep*((1.0/(xu0*xu0+xv0*xv0)) * &
194 ( (alpha*q_ts/cp_w+omg_m*beta*sss*hl_ts/le)*grav*xz0**3/(4.0*rich*rho) &
195 +( (tox*xu0+toy*xv0)/rho+(3.0*drho-alpha*i0*aw*xz0/(rho*cp_w)) &
196 *grav*xz0*xz0/(4.0*rich) )*xzts0 ))
197 xtts1 = xtts0 + timestep*(i0*aw*xzts0-q_ts)/(rho*cp_w)
198
199 ! if ( 2.0*xt1/xz1 > 0.001 ) then
200 ! write(*,'(a,i5,2f8.3,4f8.2,f10.6,10f8.4)') 'eulerm_01 : ',kdt,alat,alon,soltim/3600.,i0,q,q_warm,sep,&
201 ! 2.0*xt1/xz1,2.0*xs1/xz1,2.0*xu1/xz1,2.0*xv1/xz1,xz1,xtts1,xzts1,d_conv,t_fcl,te
202 ! endif
203
204 call sw_ps_9b(xz1,fw)
205 q_warm = fw*i0-q !total heat abs in warm layer
206 call sw_ps_9b_aw(xz1,aw)
207 drho = -alpha*q_warm/(rho*cp_w) + omg_m*beta*sep
208 dzw = xz1*(tox*xu1+toy*xv1) / (rho*(xu1*xu1+xv1*xv1)) &
209 + xz1*xz1*xz1*drho*grav / (4.0*rich*(xu1*xu1+xv1*xv1))
210
211 xt2 = xt0 + timestep*q_warm/(rho*cp_w)
212 xs2 = xs0 + timestep*sep
213 xu2 = xu0 + timestep*(fc*xv1+tox/rho)
214 xv2 = xv0 + timestep*(-fc*xu1+toy/rho)
215 xz2 = xz0 + timestep*dzw
216
217 ! if (lprnt) print *,' xt2=',xt2,' xz2=',xz2
218
219 if ( xt2 <= zero .or. xz2 <= zero .or. xz2 > z_w_max ) then
220 call dtl_reset(xt,xs,xu,xv,xz,xzts,xtts)
221 return
222 endif
223
224 xzts2 = xzts0 + timestep*((1.0/(xu1*xu1+xv1*xv1)) * &
225 ( (alpha*q_ts/cp_w+omg_m*beta*sss*hl_ts/le)*grav*xz1**3/(4.0*rich*rho) &
226 +( (tox*xu1+toy*xv1)/rho+(3.0*drho-alpha*i0*aw*xz1/(rho*cp_w))* &
227 grav*xz1*xz1/(4.0*rich) )*xzts1 ))
228 xtts2 = xtts0 + timestep*(i0*aw*xzts1-q_ts)/(rho*cp_w)
229
230 xt = 0.5*(xt1 + xt2)
231 xs = 0.5*(xs1 + xs2)
232 xu = 0.5*(xu1 + xu2)
233 xv = 0.5*(xv1 + xv2)
234 xz = 0.5*(xz1 + xz2)
235 xzts = 0.5*(xzts1 + xzts2)
236 xtts = 0.5*(xtts1 + xtts2)
237
238 if ( xt <= zero .or. xz < zero .or. xz > z_w_max ) then
239 call dtl_reset(xt,xs,xu,xv,xz,xzts,xtts)
240 endif
241
242 ! if (lprnt) print *,' xt=',xt,' xz=',xz
243 ! if ( 2.0*xt/xz > 0.001 ) then
244 ! write(*,'(a,i5,2f8.3,4f8.2,f10.6,10f8.4)') 'eulerm_02 : ',kdt,alat,alon,soltim/3600.,i0,q,q_warm,sep,&
245 ! 2.0*xt/xz,2.0*xs/xz,2.0*xu/xz,2.0*xv/xz,xz,xtts,xzts,d_conv,t_fcl,te
246 ! endif
247 return
248
249 end subroutine eulerm
250
253 subroutine dtm_1p_zwa(kdt,timestep,i0,q,rho,d_conv,xt,xs,xu,xv,xz,tr_mda,tr_fca,tr_tla,tr_mwa)
254 ! apply xz adjustment: minimum depth adjustment (mda)
255 ! free convection adjustment (fca);
256 ! top layer adjustment (tla);
257 ! maximum warming adjustment (mwa)
258 !
259 integer, intent(in) :: kdt
260 real(kind=kind_phys), intent(in) :: timestep,i0,q,rho,d_conv
261 real(kind=kind_phys), intent(inout) :: xt,xs,xu,xv,xz
262 real(kind=kind_phys), intent(out) :: tr_mda,tr_fca,tr_tla,tr_mwa
263 ! local variables
264 real(kind=kind_phys) :: dz,t0,ttop0,ttop,fw,q_warm
265 ! TODO: xz_mwa is unset but used below in max function
266 real(kind=kind_phys) :: xz_fca,xz_tla,xz_mwa
267 !
268 real(kind=kind_phys) :: xz_mda
269
270 tr_mda = zero; tr_fca = zero; tr_tla = zero; tr_mwa = zero
271
272 ! apply mda
273 if ( z_w_min > xz ) then
274 xz_mda = z_w_min
275 endif
276 ! apply fca
277 if ( d_conv > zero ) then
278 xz_fca = 2.0*xt/((2.0*xt/xz)*(1.0-d_conv/(2.0*xz)))
279 tr_fca = 1.0
280 if ( xz_fca >= z_w_max ) then
281 call dtl_reset_cv(xt,xs,xu,xv,xz)
282 go to 10
283 endif
284 endif
285 ! apply tla
286 dz = min(xz,max(d_conv,delz))
287 call sw_ps_9b(dz,fw)
288 q_warm=fw*i0-q !total heat abs in warm layer
289
290 if ( q_warm > zero ) then
291 call cal_ttop(kdt,timestep,q_warm,rho,dz,xt,xz,ttop0)
292 ! ttop = (2.0*xt/xz)*(1.0-dz/(2.0*xz))
293 ttop = ((xt+xt)/xz)*(1.0-dz/(xz+xz))
294 if ( ttop > ttop0 ) then
295 xz_tla = (xt+sqrt(xt*(xt-delz*ttop0)))/ttop0
296 tr_tla = 1.0
297 if ( xz_tla >= z_w_max ) then
298 call dtl_reset_cv(xt,xs,xu,xv,xz)
299 go to 10
300 endif
301 endif
302 endif
303
304 ! apply mwa
305 t0 = 2.0*xt/xz
306 if ( t0 > tw_max ) then
307 if ( xz >= z_w_max ) then
308 call dtl_reset_cv(xt,xs,xu,xv,xz)
309 go to 10
310 endif
311 endif
312
313 xz = max(xz_mda,xz_fca,xz_tla,xz_mwa)
314
31510 continue
316
317 end subroutine dtm_1p_zwa
318
321 subroutine dtm_1p_fca(d_conv,xt,xtts,xz,xzts)
322
323 ! apply xz adjustment: free convection adjustment (fca);
324 !
325 real(kind=kind_phys), intent(in) :: d_conv,xt,xtts
326 real(kind=kind_phys), intent(inout) :: xz,xzts
327 ! local variables
328 real(kind=kind_phys) :: t_fcl,t0
329 !
330 t0 = 2.0*xt/xz
331 t_fcl = t0*(1.0-d_conv/(2.0*xz))
332 xz = 2.0*xt/t_fcl
333 ! xzts = 2.0*xtts/t_fcl
334
335 end subroutine dtm_1p_fca
336
339 subroutine dtm_1p_tla(dz,te,xt,xtts,xz,xzts)
340
341 ! apply xz adjustment: top layer adjustment (tla);
342 !
343 real(kind=kind_phys), intent(in) :: dz,te,xt,xtts
344 real(kind=kind_phys), intent(inout) :: xz,xzts
345 ! local variables
346 real(kind=kind_phys) :: tem
347 !
348 tem = xt*(xt-dz*te)
349 if (tem > zero) then
350 xz = (xt+sqrt(xt*(xt-dz*te)))/te
351 else
352 xz = z_w_max
353 endif
354 ! xzts = xtts*(1.0+0.5*(2.0*xt-dz*te)/sqrt(xt*(xt-dz*te)))/te
355 end subroutine dtm_1p_tla
356
359 subroutine dtm_1p_mwa(xt,xtts,xz,xzts)
360
361 ! apply xz adjustment: maximum warming adjustment (mwa)
362 !
363 real(kind=kind_phys), intent(in) :: xt,xtts
364 real(kind=kind_phys), intent(inout) :: xz,xzts
365 ! local variables
366 !
367 xz = 2.0*xt/tw_max
368 ! xzts = 2.0*xtts/tw_max
369 end subroutine dtm_1p_mwa
370
373 subroutine dtm_1p_mda(xt,xtts,xz,xzts)
374
375 ! apply xz adjustment: minimum depth adjustment (mda)
376 !
377 real(kind=kind_phys), intent(in) :: xt,xtts
378 real(kind=kind_phys), intent(inout) :: xz,xzts
379 ! local variables
380 real(kind=kind_phys) :: ta
381 !
382 xz = max(z_w_min,xz)
383 ta = 2.0*xt/xz
384 ! xzts = 2.0*xtts/ta
385
386 end subroutine dtm_1p_mda
387
390 subroutine dtm_1p_mta(dta,xt,xtts,xz,xzts)
391
392 ! apply xz adjustment: maximum temperature adjustment (mta)
393 !
394 real(kind=kind_phys), intent(in) :: dta,xt,xtts
395 real(kind=kind_phys), intent(inout) :: xz,xzts
396 ! local variables
397 real(kind=kind_phys) :: ta
398 !
399 ta = max(zero,2.0*xt/xz-dta)
400 if ( ta > zero ) then
401 xz = 2.0*xt/ta
402 else
403 xz = z_w_max
404 endif
405 ! xzts = 2.0*xtts/ta
406
407 end subroutine dtm_1p_mta
408
411 subroutine convdepth(kdt,timestep,i0,q,sss,sep,rho,alpha,beta,xt,xs,xz,d_conv)
412
413 !
414 ! calculate depth for convective adjustment
415 !
416
417 integer, intent(in) :: kdt
418 real(kind=kind_phys), intent(in) :: timestep,i0,q,sss,sep,rho,alpha,beta
419 real(kind=kind_phys), intent(in) :: xt,xs,xz
420 real(kind=kind_phys), intent(out) :: d_conv
421 real(kind=kind_phys) :: t,s,d_conv_ini,d_conv2,fxp,aw,s1,s2,fac1
422 integer :: n
423 !
424 ! input variables
425 !
426 ! timestep: time step in seconds
427 ! i0 : solar radiation flux at the surface (wm^-2)
428 ! q : non-solar heat flux at the surface (wm^-2)
429 ! sss : salinity (ppt)
430 ! sep : sr(e-p) (ppt*m/s)
431 ! rho : sea water density (kg*m^-3)
432 ! alpha : thermal expansion coefficient (1/k)
433 ! beta : saline contraction coefficient (1/ppt)
434 ! xt : initial heat content (k*m)
435 ! xs : initial salinity content (ppt*m)
436 ! xz : initial dtl thickness (m)
437 !
438 ! output variables
439 !
440 ! d_conv : free convection depth (m)
441
442 ! t : initial diurnal warming t (k)
443 ! s : initial diurnal warming s (ppt)
444
445 n = 0
446 t = 2.0*xt/xz
447 s = 2.0*xs/xz
448
449 s1 = alpha*rho*t-omg_m*beta*rho*s
450
451 if ( s1 == zero ) then
452 d_conv = zero
453 else
454
455 fac1 = alpha*q/cp_w+omg_m*beta*rho*sep
456 if ( i0 <= zero ) then
457 d_conv2=(2.0*xz*timestep/s1)*fac1
458 if ( d_conv2 > zero ) then
459 d_conv = sqrt(d_conv2)
460 else
461 d_conv = zero
462 endif
463 elseif ( i0 > zero ) then
464
465 d_conv_ini = zero
466
467 iter_conv: do n = 1, niter_conv
468 call sw_ps_9b(d_conv_ini,fxp)
469 call sw_ps_9b_aw(d_conv_ini,aw)
470 s2 = alpha*(q-(fxp-aw*d_conv_ini)*i0)/cp_w+omg_m*beta*rho*sep
471 d_conv2=(2.0*xz*timestep/s1)*s2
472 if ( d_conv2 < zero ) then
473 d_conv = zero
474 exit iter_conv
475 endif
476 d_conv = sqrt(d_conv2)
477 if ( abs(d_conv-d_conv_ini) < eps_conv .and. n <= niter_conv ) exit iter_conv
478 d_conv_ini = d_conv
479 enddo iter_conv
480 d_conv = max(zero,min(d_conv,z_w_max))
481 endif ! if ( i0 <= zero ) then
482
483 endif ! if ( s1 == zero ) then
484
485 ! if ( d_conv > 0.01 ) then
486 ! write(*,'(a,i4,i3,10f9.3,3f10.6,f10.1,f6.2)') ' d_conv : ',kdt,n,d_conv,d_conv_ini,q,i0,rho,cp_w,timestep,xt,xs,xz,sep, &
487 ! s1,s2,d_conv2,aw
488 ! endif
489
490 end subroutine convdepth
491
493 subroutine dtm_onset(kdt,timestep,rich,tox,toy,i0,q,sss,sep,q_ts,hl_ts,rho, &
494 alpha,beta,alon,sinlat,soltim,grav,le,xt,xs,xu,xv,xz,xzts,xtts)
495 !
496 ! determine xz iteratively (starting wit fw = 0.5) and then update the other 6 variables
497 !
498
499 integer,intent(in) :: kdt
500 real(kind=kind_phys), intent(in) :: timestep,rich,tox,toy,i0,q,sss,sep,q_ts, &
501 hl_ts,rho,alpha,beta,alon,sinlat,soltim,grav,le
502 real(kind=kind_phys), intent(out) :: xt,xs,xu,xv,xz,xzts,xtts
503 real(kind=kind_phys) :: xt0,xs0,xu0,xv0,xz0
504 real(kind=kind_phys) :: xt1,xs1,xu1,xv1,xz1
505 real(kind=kind_phys) :: fw,aw,q_warm,ft0,fs0,fu0,fv0,fz0,ft1,fs1,fu1,fv1,fz1
506 real(kind=kind_phys) :: coeff1,coeff2,ftime,z_w,z_w_tmp,fc,warml,alat
507 integer :: n
508 !
509 ! input variables
510 !
511 ! timestep: time step in seconds
512 ! tox : x wind stress (n*m^-2 or kg/m/s^2)
513 ! toy : y wind stress (n*m^-2 or kg/m/s^2)
514 ! i0 : solar radiation flux at the surface (wm^-2)
515 ! q : non-solar heat flux at the surface (wm^-2)
516 ! sss : salinity (ppt)
517 ! sep : sr(e-p) (ppt*m/s)
518 ! rho : sea water density (kg*m^-3)
519 ! alpha : thermal expansion coefficient (1/k)
520 ! beta : saline contraction coefficient (1/ppt)
521 ! alon : longitude
522 ! sinlat : sine(latitude)
523 ! grav : gravity accelleration
524 ! le : le=(2.501-.00237*tsea)*1e6
525 !
526 ! output variables
527 !
528 ! xt : onset t content in dtl
529 ! xs : onset s content in dtl
530 ! xu : onset u content in dtl
531 ! xv : onset v content in dtl
532 ! xz : onset dtl thickness (m)
533 ! xzts : onset d(xz)/d(ts) (m/k )
534 ! xtts : onset d(xt)/d(ts) (m)
535
536 fc=1.46/10000.0/2.0*sinlat
537 alat = asin(sinlat)
538 !
539 ! initializing dtl (just before the onset)
540 !
541 xt0 = zero
542 xs0 = zero
543 xu0 = zero
544 xv0 = zero
545
546 z_w_tmp=z_w_ini
547
548 call sw_ps_9b(z_w_tmp,fw)
549 ! fw=0.5 !
550 q_warm=fw*i0-q !total heat abs in warm layer
551
552 if ( abs(alat) > 1.0 ) then
553 ftime=sqrt((2.0-2.0*cos(fc*timestep))/(fc*fc*timestep))
554 else
555 ftime=timestep
556 endif
557
558 coeff1=alpha*grav/cp_w
559 coeff2=omg_m*beta*grav*rho
560 warml = coeff1*q_warm-coeff2*sep
561
562 if ( warml > zero .and. q_warm > zero) then
563 iters_z_w: do n = 1,niter_z_w
564 if ( warml > zero .and. q_warm > zero ) then
565 z_w=sqrt(2.0*rich*ftime/rho)*sqrt(tox**2+toy**2)/sqrt(warml)
566 else
567 z_w = z_w_max
568 exit iters_z_w
569 endif
570
571 ! write(*,'(a,i4,i4,10f9.3,f9.6,f3.0)') ' z_w = ',kdt,n,z_w,z_w_tmp,timestep,q_warm,q,i0,fw,tox,toy,sep,warml,omg_m
572
573 if (abs(z_w - z_w_tmp) < eps_z_w .and. z_w/=z_w_max .and. n < niter_z_w) exit iters_z_w
574 z_w_tmp=z_w
575 call sw_ps_9b(z_w_tmp,fw)
576 q_warm = fw*i0-q
577 warml = coeff1*q_warm-coeff2*sep
578 end do iters_z_w
579 else
580 z_w=z_w_max
581 endif
582
583 xz0 = max(z_w,z_w_min)
584
585 !
586 ! update xt, xs, xu, xv
587 !
588 if ( z_w < z_w_max .and. q_warm > zero) then
589
590 call sw_ps_9b(z_w,fw)
591 q_warm=fw*i0-q !total heat abs in warm layer
592
593 ft0 = q_warm/(rho*cp_w)
594 fs0 = sep
595 fu0 = fc*xv0+tox/rho
596 fv0 = -fc*xu0+toy/rho
597
598 xt1 = xt0 + timestep*ft0
599 xs1 = xs0 + timestep*fs0
600 xu1 = xu0 + timestep*fu0
601 xv1 = xv0 + timestep*fv0
602
603 fz0 = xz0*((tox*xu1+toy*xv1)/rho+omg_m*beta*grav*sep*xz0*xz0/(4.0*rich) &
604 -alpha*grav*q_warm*xz0*xz0/(4.0*rich*cp_w*rho))/(xu1*xu1+xv1*xv1)
605 xz1 = xz0 + timestep*fz0
606
607 xz1 = max(xz1,z_w_min)
608
609 if ( xt1 < zero .or. xz1 > z_w_max ) then
610 call dtl_reset(xt,xs,xu,xv,xz,xtts,xzts)
611 return
612 endif
613
614 call sw_ps_9b(xz1,fw)
615 q_warm=fw*i0-q !total heat abs in warm layer
616
617 ft1 = q_warm/(rho*cp_w)
618 fs1 = sep
619 fu1 = fc*xv1+tox/rho
620 fv1 = -fc*xu1+toy/rho
621
622 fz1 = xz1*((tox*xu1+toy*xv1)/rho+omg_m*beta*grav*sep*xz1*xz1/(4.0*rich) &
623 -alpha*grav*q_warm*xz1*xz1/(4.0*rich*cp_w*rho))/(xu1*xu1+xv1*xv1)
624
625 xt = xt0 + 0.5*timestep*(ft0+ft1)
626 xs = xs0 + 0.5*timestep*(fs0+fs1)
627 xu = xu0 + 0.5*timestep*(fu0+fu1)
628 xv = xv0 + 0.5*timestep*(fv0+fv1)
629 xz = xz0 + 0.5*timestep*(fz0+fz1)
630
631 xz = max(xz,z_w_min)
632
633 call sw_ps_9b_aw(xz,aw)
634
635 ! xzts = (q_ts+(cp_w*omg_m*beta*sss/(le*alpha))*hl_ts)*xz/(i0*xz*aw+2.0*q_warm-2.0*(rho*cp_w*omg_m*beta*sss/alpha)*(sep/sss))
636 xzts = (q_ts+omg_m*rho*cp_w*beta*sss*hl_ts*xz/(le*alpha))/(i0*xz*aw+2.0*q_warm-2.0*omg_m*rho*cp_w*beta*sss*sep/(le*alpha))
637 xtts = timestep*(i0*aw*xzts-q_ts)/(rho*cp_w)
638 endif
639
640 if ( xt < zero .or. xz > z_w_max ) then
641 call dtl_reset(xt,xs,xu,xv,xz,xtts,xzts)
642 endif
643
644 return
645
646 end subroutine dtm_onset
647
651 subroutine cal_w(kdt,xz,xt,xzts,xtts,w_0,w_d)
652 !
653 ! abstract: calculate w_0,w_d
654 !
655 ! input variables
656 !
657 ! kdt : the number of time step
658 ! xt : dtl heat content
659 ! xz : dtl depth
660 ! xzts : d(zw)/d(ts)
661 ! xtts : d(xt)/d(ts)
662 !
663 ! output variables
664 !
665 ! w_0 : coefficint 1 to calculate d(tw)/d(ts)
666 ! w_d : coefficint 2 to calculate d(tw)/d(ts)
667
668 integer, intent(in) :: kdt
669 real(kind=kind_phys), intent(in) :: xz,xt,xzts,xtts
670 real(kind=kind_phys), intent(out) :: w_0,w_d
671
672 w_0 = 2.0*(xtts-xt*xzts/xz)/xz
673 w_d = (2.0*xt*xzts/xz**2-w_0)/xz
674
675 ! if ( 2.0*xt/xz > 1.0 ) then
676 ! write(*,'(a,i4,2f9.3,4f10.4))') ' cal_w : ',kdt,xz,xt,w_0,w_d,xzts,xtts
677 ! endif
678 end subroutine cal_w
679
683 subroutine cal_ttop(kdt,timestep,q_warm,rho,dz,xt,xz,ttop)
684 !
685 ! abstract: calculate
686 !
687 ! input variables
688 !
689 ! kdt : the number of record
690 ! timestep : the number of record
691 ! q_warm : total heat abs in layer dz
692 ! rho : sea water density
693 ! dz : dz = max(delz,d_conv) top layer thickness defined to adjust xz
694 ! xt : heat content in dtl at previous time
695 ! xz : dtl thickness at previous time
696 !
697 ! output variables
698 !
699 ! ttop : the diurnal warming amount at the top layer with thickness of delz
700
701 integer, intent(in) :: kdt
702 real(kind=kind_phys), intent(in) :: timestep,q_warm,rho,dz,xt,xz
703 real(kind=kind_phys), intent(out) :: ttop
704 real(kind=kind_phys) :: dt_warm,t0
705
706 dt_warm = (xt+xt)/xz
707 t0 = dt_warm*(1.0-dz/(xz+xz))
708 ttop = t0 + q_warm*timestep/(rho*cp_w*dz)
709
710 end subroutine cal_ttop
711
715 subroutine app_sfs(kdt,xt,xs,xu,xv,alpha,beta,grav,d_1p,xz)
716 !
717 ! abstract: adjust dtm-1p dtl thickness by applying shear flow stability with assumed exponetial profile
718 !
719 ! input variables
720 !
721 ! kdt : the number of record
722 ! xt : heat content in dtl
723 ! xs : salinity content in dtl
724 ! xu : u-current content in dtl
725 ! xv : v-current content in dtl
726 ! alpha
727 ! beta
728 ! grav
729 ! d_1p : dtl depth before sfs applied
730 !
731 ! output variables
732 !
733 ! xz : dtl depth
734
735 integer, intent(in) :: kdt
736 real(kind=kind_phys), intent(in) :: xt,xs,xu,xv,alpha,beta,grav,d_1p
737 real(kind=kind_phys), intent(out) :: xz
738 ! real(kind=kind_phys) :: ze,cc,xz0,l,d_sfs, t_sfs, tem
739 real(kind=kind_phys) :: cc,l,d_sfs,tem
740 real(kind=kind_phys), parameter :: c2 = 0.3782
741
742 cc = ri_g/(grav*c2)
743
744 tem = alpha*xt - beta*xs
745 if (tem > zero) then
746 d_sfs = sqrt(2.0*cc*(xu*xu+xv*xv)/tem)
747 else
748 d_sfs = zero
749 endif
750
751 ! xz0 = d_1p
752 ! iter_sfs: do n = 1, niter_sfs
753 ! l = int_epn(0.0,xz0,0.0,xz0,2)
754 ! d_sfs = cc*(xu*xu+xv*xv)/((alpha*xt-beta*xs)*l)
755 ! write(*,'(a,i6,i3,4f9.4))') ' app_sfs_iter : ',kdt,n,cc,l,xz0,d_sfs
756 ! if ( abs(d_sfs-xz0) < eps_sfs .and. n <= niter_sfs ) exit iter_sfs
757 ! xz0 = d_sfs
758 ! enddo iter_sfs
759
760 ! ze = a2*d_sfs ! not used!
761
762 l = int_epn(zero,d_sfs,zero,d_sfs,2)
763
764 ! t_sfs = xt/l
765 ! xz = (xt+xt) / t_sfs
766
767 xz = l + l
768
769 ! write(*,'(a,i6,6f9.4))') ' app_sfs : ',kdt,xz0,d_sfs,d_1p,xz,2.0*xt/d_1p,t_sfs
770 end subroutine app_sfs
771
774 subroutine cal_tztr(kdt,xt,c_0,c_d,w_0,w_d,zc,zw,z,tztr)
775 !
776 ! abstract: calculate d(tz)/d(ts)
777 !
778 ! input variables
779 !
780 ! kdt : the number of record
781 ! xt : heat content in dtl
782 ! xz : dtl depth (m)
783 ! c_0 : coefficint 1 to calculate d(tc)/d(ts)
784 ! c_d : coefficint 2 to calculate d(tc)/d(ts)
785 ! w_0 : coefficint 1 to calculate d(tw)/d(ts)
786 ! w_d : coefficint 2 to calculate d(tw)/d(ts)
787 !
788 ! output variables
789 !
790 ! tztr : d(tz)/d(tr)
791
792 integer, intent(in) :: kdt
793 real(kind=kind_phys), intent(in) :: xt,c_0,c_d,w_0,w_d,zc,zw,z
794 real(kind=kind_phys), intent(out) :: tztr
795
796 if ( xt > zero ) then
797 if ( z <= zc ) then
798 ! tztr = 1.0/(1.0-w_0+c_0)+z*(w_d-c_d)/(1.0-w_0+c_0)
799 tztr = (1.0+z*(w_d-c_d))/(1.0-w_0+c_0)
800 elseif ( z > zc .and. z < zw ) then
801 ! tztr = (1.0+c_0)/(1.0-w_0+c_0)+z*w_d/(1.0-w_0+c_0)
802 tztr = (1.0+c_0+z*w_d)/(1.0-w_0+c_0)
803 elseif ( z >= zw ) then
804 tztr = 1.0
805 endif
806 elseif ( xt == zero ) then
807 if ( z <= zc ) then
808 ! tztr = 1.0/(1.0+c_0)-z*c_d/(1.0+c_0)
809 tztr = (1.0-z*c_d)/(1.0+c_0)
810 else
811 tztr = 1.0
812 endif
813 else
814 tztr = 1.0
815 endif
816
817 ! write(*,'(a,i4,9f9.4))') ' cal_tztr : ',kdt,xt,c_0,c_d,w_0,w_d,zc,zw,z,tztr
818 end subroutine cal_tztr
819
823 subroutine 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)
824 !
825 ! upper ocean cool-skin parameterizaion, fairall et al, 1996.
826 !
827 ! input:
828 ! ustar_a : atmosphreic friction velocity at the air-sea interface (m/s)
829 ! f_nsol : the "nonsolar" part of the surface heat flux (w/m^s)
830 ! f_sol_0 : solar radiation at the ocean surface (w/m^2)
831 ! evap : latent heat flux (w/m^2)
832 ! sss : ocean upper mixed layer salinity (ppu)
833 ! alpha : thermal expansion coefficient
834 ! beta : saline contraction coefficient
835 ! rho_w : oceanic density
836 ! rho_a : atmospheric density
837 ! ts : oceanic surface temperature
838 ! q_ts : d(q)/d(ts) : q = the sum of non-solar heat fluxes
839 ! hl_ts : d(hl)/d(ts)
840 ! grav : gravity
841 ! le :
842 !
843 ! output:
844 ! deltat_c: cool-skin temperature correction (degrees k)
845 ! z_c : molecular sublayer (cool-skin) thickness (m)
846 ! c_0 : coefficient1 to calculate d(tz)/d(ts)
847 ! c_d : coefficient2 to calculate d(tz)/d(ts)
848
849 !
850 real(kind=kind_phys), intent(in) :: ustar_a,f_nsol,f_sol_0,evap,sss,alpha,beta,rho_w,rho_a,ts,q_ts,hl_ts,grav,le
851 real(kind=kind_phys), intent(out) :: deltat_c,z_c,c_0,c_d
852 ! declare local variables
853 real(kind=kind_phys), parameter :: a1=0.065, a2=11.0, a3=6.6e-5, a4=8.0e-4, tcw=0.6, tcwi=1.0/tcw
854 real(kind=kind_phys) :: a_c,b_c,zc_ts,bc1,bc2
855 real(kind=kind_phys) :: xi,hb,ustar1_a,bigc,deltaf,fxp
856 real(kind=kind_phys) :: zcsq
857 real(kind=kind_phys) :: cc1,cc2,cc3
858
859
860 z_c = z_c_ini ! initial guess
861
862 ustar1_a = max(ustar_a,ustar_a_min)
863
864 call sw_rad_skin(z_c,fxp)
865 deltaf = f_sol_0*fxp
866
867 hb = alpha*(f_nsol-deltaf)+beta*sss*cp_w*evap/le
868 bigc = 16*grav*cp_w*(rho_w*visw)**3/(rho_a*rho_a*kw*kw)
869
870 if ( hb > 0 ) then
871 xi = 6./(1+(bigc*hb/ustar1_a**4)**0.75)**0.3333333
872 else
873 xi = 6.0
874 endif
875 z_c = min(z_c_max,xi*visw/(sqrt(rho_a/rho_w)*ustar1_a ))
876
877 call sw_rad_skin(z_c,fxp)
878
879 deltaf = f_sol_0*fxp
880 deltaf = f_nsol - deltaf
881 if ( deltaf > 0 ) then
882 deltat_c = deltaf * z_c / kw
883 else
884 deltat_c = zero
885 z_c = zero
886 endif
887 !
888 ! calculate c_0 & c_d
889 !
890 if ( z_c > zero ) then
891 cc1 = 6.0*visw / (tcw*ustar1_a*sqrt(rho_a/rho_w))
892 cc2 = bigc*alpha / max(ustar_a,ustar_a_min)**4
893 cc3 = beta*sss*cp_w/(alpha*le)
894 zcsq = z_c * z_c
895 a_c = a2 + a3/zcsq - (a3/(a4*z_c)+a3/zcsq) * exp(-z_c/a4)
896
897 if ( hb > zero .and. zcsq > zero .and. alpha > zero) then
898 bc1 = zcsq * (q_ts+cc3*hl_ts)
899 bc2 = zcsq * f_sol_0*a_c - 4.0*(cc1*tcw)**3*(hb/alpha)**0.25/(cc2**0.75*zcsq)
900 zc_ts = bc1/bc2
901 ! b_c = z_c**2*(q_ts+cc3*hl_ts)/(z_c**2*f_sol_0*a_c-4.0*(cc1*tcw)**3*(hb/alpha)**0.25/(cc2**0.75*z_c**2)) ! d(z_c)/d(ts)
902 b_c = (q_ts+cc3*hl_ts)/(f_sol_0*a_c &
903 - 4.0*(cc1*tcw)**3*(hb/alpha)**0.25/(cc2**0.75*zcsq*zcsq)) ! d(z_c)/d(ts)
904 c_0 = (z_c*q_ts+(f_nsol-deltaf-f_sol_0*a_c*z_c)*b_c)*tcwi
905 c_d = (f_sol_0*a_c*z_c*b_c-q_ts)*tcwi
906
907 else
908 b_c = zero
909 zc_ts = zero
910 c_0 = z_c*q_ts*tcwi
911 c_d = -q_ts*tcwi
912 endif
913
914 ! if ( c_0 < 0.0 ) then
915 ! write(*,'(a,2f12.6,10f10.6)') ' c_0, c_d = ',c_0,c_d,b_c,zc_ts,hb,bc1,bc2,z_c,cc1,cc2,cc3,z_c**2
916 ! endif
917
918 ! c_0 = z_c*q_ts/tcw
919 ! c_d = -q_ts/tcw
920
921 else
922 c_0 = zero
923 c_d = zero
924 endif ! if ( z_c > 0.0 ) then
925
926 end subroutine cool_skin
927 !
928 !======================
929 !
932 real function int_epn(z1,z2,zmx,ztr,n)
933 !
934 ! abstract: calculate a definitive integral of an exponetial curve (power of 2)
935 !
936 real(kind_phys) :: z1,z2,zmx,ztr,zi
937 real(kind_phys) :: fa,fb,fi,int
938 integer :: m,i,n
939
940 m = nint((z2-z1)/delz)
941 fa = exp(-exp_const*((z1-zmx)/(ztr-zmx))**n)
942 fb = exp(-exp_const*((z2-zmx)/(ztr-zmx))**n)
943 int = zero
944 do i = 1, m-1
945 zi = z1 + delz*float(i)
946 fi = exp(-exp_const*((zi-zmx)/(ztr-zmx))**n)
947 int = int + fi
948 enddo
949 int_epn = delz*((fa+fb)/2.0 + int)
950 end function int_epn
951
954 subroutine dtl_reset_cv(xt,xs,xu,xv,xz)
955 real(kind=kind_phys), intent(inout) :: xt,xs,xu,xv,xz
956 xt = zero
957 xs = zero
958 xu = zero
959 xv = zero
960 xz = z_w_max
961 end subroutine dtl_reset_cv
962
965 subroutine dtl_reset(xt,xs,xu,xv,xz,xzts,xtts)
966 real(kind=kind_phys), intent(inout) :: xt,xs,xu,xv,xz,xzts,xtts
967 xt = zero
968 xs = zero
969 xu = zero
970 xv = zero
971 xz = z_w_max
972 xtts = zero
973 xzts = zero
974 end subroutine dtl_reset
975
976end module nst_module
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 dtm_1p_zwa(kdt, timestep, i0, q, rho, d_conv, xt, xs, xu, xv, xz, tr_mda, tr_fca, tr_tla, tr_mwa)
This subroutine applies xz adjustment.
subroutine, public dtm_1p_fca(d_conv, xt, xtts, xz, xzts)
This subroutine applies free convection adjustment(fca).
real function int_epn(z1, z2, zmx, ztr, n)
This function calculates a definitive integral of an exponential curve (power of 2).
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 app_sfs(kdt, xt, xs, xu, xv, alpha, beta, grav, d_1p, xz)
This subroutine adjust dtm-1p dtl thickness by applying shear flow stability with assumed exponential...
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 dtm_onset(kdt, timestep, rich, tox, toy, i0, q, sss, sep, q_ts, hl_ts, rho, alpha, beta, alon, sinlat, soltim, grav, le, xt, xs, xu, xv, xz, xzts, xtts)
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 eulerm(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 integrates one time step with modified Euler method.
subroutine dtl_reset_cv(xt, xs, xu, xv, xz)
This subroutine resets the value of xt,xs,xu,xv,xz.
subroutine cal_tztr(kdt, xt, c_0, c_d, w_0, w_d, zc, zw, z, tztr)
This subroutine calculates d(tz)/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,...
This module contains the diurnal thermocline layer model (DTM) of the GFS NSST scheme.