16 use machine ,
only : kind_phys
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)
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,&
44 real(kind=kind_phys),
intent(inout) :: xt,xs,xu,xv,xz,xzts,xtts
81 if ( xt <= zero )
then
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
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)
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)
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, &
108 real(kind=kind_phys),
intent(inout) :: xt,xs,xu,xv,xz,xzts,xtts
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
158 speed = max(1.0e-8, xu0*xu0+xv0*xv0)
160 alat = asin(sinlat)*rad2deg
162 fc = const_rot*sinlat
170 drho = -alpha*q_warm/(rho*cp_w) + omg_m*beta*sep
174 dzw = xz0 * ((tox*xu0+toy*xv0) / (rho*speed) &
175 + xz0*xz0*drho*grav / (4.0*rich*speed))
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
186 if ( xt1 <= zero .or. xz1 <= zero .or. xz1 > z_w_max )
then
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)
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))
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
219 if ( xt2 <= zero .or. xz2 <= zero .or. xz2 > z_w_max )
then
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)
235 xzts = 0.5*(xzts1 + xzts2)
236 xtts = 0.5*(xtts1 + xtts2)
238 if ( xt <= zero .or. xz < zero .or. xz > z_w_max )
then
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)
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
264 real(kind=kind_phys) :: dz,t0,ttop0,ttop,fw,q_warm
266 real(kind=kind_phys) :: xz_fca,xz_tla,xz_mwa
268 real(kind=kind_phys) :: xz_mda
270 tr_mda = zero; tr_fca = zero; tr_tla = zero; tr_mwa = zero
273 if ( z_w_min > xz )
then
277 if ( d_conv > zero )
then
278 xz_fca = 2.0*xt/((2.0*xt/xz)*(1.0-d_conv/(2.0*xz)))
280 if ( xz_fca >= z_w_max )
then
286 dz = min(xz,max(d_conv,delz))
290 if ( q_warm > zero )
then
291 call cal_ttop(kdt,timestep,q_warm,rho,dz,xt,xz,ttop0)
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
297 if ( xz_tla >= z_w_max )
then
306 if ( t0 > tw_max )
then
307 if ( xz >= z_w_max )
then
313 xz = max(xz_mda,xz_fca,xz_tla,xz_mwa)
325 real(kind=kind_phys),
intent(in) :: d_conv,xt,xtts
326 real(kind=kind_phys),
intent(inout) :: xz,xzts
328 real(kind=kind_phys) :: t_fcl,t0
331 t_fcl = t0*(1.0-d_conv/(2.0*xz))
343 real(kind=kind_phys),
intent(in) :: dz,te,xt,xtts
344 real(kind=kind_phys),
intent(inout) :: xz,xzts
346 real(kind=kind_phys) :: tem
350 xz = (xt+sqrt(xt*(xt-dz*te)))/te
363 real(kind=kind_phys),
intent(in) :: xt,xtts
364 real(kind=kind_phys),
intent(inout) :: xz,xzts
377 real(kind=kind_phys),
intent(in) :: xt,xtts
378 real(kind=kind_phys),
intent(inout) :: xz,xzts
380 real(kind=kind_phys) :: ta
394 real(kind=kind_phys),
intent(in) :: dta,xt,xtts
395 real(kind=kind_phys),
intent(inout) :: xz,xzts
397 real(kind=kind_phys) :: ta
399 ta = max(zero,2.0*xt/xz-dta)
400 if ( ta > zero )
then
411 subroutine convdepth(kdt,timestep,i0,q,sss,sep,rho,alpha,beta,xt,xs,xz,d_conv)
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
449 s1 = alpha*rho*t-omg_m*beta*rho*s
451 if ( s1 == zero )
then
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)
463 elseif ( i0 > zero )
then
467 iter_conv:
do n = 1, niter_conv
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
476 d_conv = sqrt(d_conv2)
477 if ( abs(d_conv-d_conv_ini) < eps_conv .and. n <= niter_conv )
exit iter_conv
480 d_conv = max(zero,min(d_conv,z_w_max))
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)
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
536 fc=1.46/10000.0/2.0*sinlat
552 if ( abs(alat) > 1.0 )
then
553 ftime=sqrt((2.0-2.0*cos(fc*timestep))/(fc*fc*timestep))
558 coeff1=alpha*grav/cp_w
559 coeff2=omg_m*beta*grav*rho
560 warml = coeff1*q_warm-coeff2*sep
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)
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
577 warml = coeff1*q_warm-coeff2*sep
583 xz0 = max(z_w,z_w_min)
588 if ( z_w < z_w_max .and. q_warm > zero)
then
593 ft0 = q_warm/(rho*cp_w)
596 fv0 = -fc*xu0+toy/rho
598 xt1 = xt0 + timestep*ft0
599 xs1 = xs0 + timestep*fs0
600 xu1 = xu0 + timestep*fu0
601 xv1 = xv0 + timestep*fv0
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
607 xz1 = max(xz1,z_w_min)
609 if ( xt1 < zero .or. xz1 > z_w_max )
then
617 ft1 = q_warm/(rho*cp_w)
620 fv1 = -fc*xu1+toy/rho
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)
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)
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)
640 if ( xt < zero .or. xz > z_w_max )
then
651 subroutine cal_w(kdt,xz,xt,xzts,xtts,w_0,w_d)
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
672 w_0 = 2.0*(xtts-xt*xzts/xz)/xz
673 w_d = (2.0*xt*xzts/xz**2-w_0)/xz
683 subroutine cal_ttop(kdt,timestep,q_warm,rho,dz,xt,xz,ttop)
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
707 t0 = dt_warm*(1.0-dz/(xz+xz))
708 ttop = t0 + q_warm*timestep/(rho*cp_w*dz)
715 subroutine app_sfs(kdt,xt,xs,xu,xv,alpha,beta,grav,d_1p,xz)
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
739 real(kind=kind_phys) :: cc,l,d_sfs,tem
740 real(kind=kind_phys),
parameter :: c2 = 0.3782
744 tem = alpha*xt - beta*xs
746 d_sfs = sqrt(2.0*cc*(xu*xu+xv*xv)/tem)
762 l =
int_epn(zero,d_sfs,zero,d_sfs,2)
774 subroutine cal_tztr(kdt,xt,c_0,c_d,w_0,w_d,zc,zw,z,tztr)
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
796 if ( xt > zero )
then
799 tztr = (1.0+z*(w_d-c_d))/(1.0-w_0+c_0)
800 elseif ( z > zc .and. z < zw )
then
802 tztr = (1.0+c_0+z*w_d)/(1.0-w_0+c_0)
803 elseif ( z >= zw )
then
806 elseif ( xt == zero )
then
809 tztr = (1.0-z*c_d)/(1.0+c_0)
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)
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
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
862 ustar1_a = max(ustar_a,ustar_a_min)
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)
871 xi = 6./(1+(bigc*hb/ustar1_a**4)**0.75)**0.3333333
875 z_c = min(z_c_max,xi*visw/(sqrt(rho_a/rho_w)*ustar1_a ))
880 deltaf = f_nsol - deltaf
881 if ( deltaf > 0 )
then
882 deltat_c = deltaf * z_c / kw
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)
895 a_c = a2 + a3/zcsq - (a3/(a4*z_c)+a3/zcsq) * exp(-z_c/a4)
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)
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))
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
936 real(kind_phys) :: z1,z2,zmx,ztr,zi
937 real(kind_phys) :: fa,fb,fi,int
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)
945 zi = z1 + delz*float(i)
946 fi = exp(-exp_const*((zi-zmx)/(ztr-zmx))**n)
949 int_epn = delz*((fa+fb)/2.0 + int)
955 real(kind=kind_phys),
intent(inout) :: xt,xs,xu,xv,xz
966 real(kind=kind_phys),
intent(inout) :: xt,xs,xu,xv,xz,xzts,xtts
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.