17 subroutine cires_ugwpv1_ngw_solv2(mpi_id, master, im, levs, kdt, dtp, &
18 tau_ngw, tm , um, vm, qm, prsl, prsi, zmet, zmeti, prslk, &
19 xlatd, sinlat, coslat, &
20 pdudt, pdvdt, pdtdt, dked, zngw)
30 use machine,
only : kind_phys
40 use ugwp_common ,
only : rgrav, grav, cpd, rd, rv, rcpdl, grav2cpd, &
41 omega2, rcpd, rcpd2, pi, pi2, fv, &
42 rad_to_deg, deg_to_rad, &
43 rdi, gor, grcp, gocp, &
44 bnv2min, bnv2max, dw2min, velmin, gr2, &
45 hpscale, rhp, rh4, grav2, rgrav2, mkzmin, mkz2min
47 use ugwp_wmsdis_init,
only : v_kxw, rv_kxw, v_kxw2, tamp_mpa, tau_min, ucrit, &
49 nslope, ilaunch, zms, &
50 zci, zdci, zci4, zci3, zci2, &
51 zaz_fct, zcosang, zsinang, nwav, nazd, &
52 zcimin, zcimax, rimin, sc2, sc2u, ric
56 real(kind=kind_phys),
parameter :: zsp_gw = 106.5e3
57 real(kind=kind_phys),
parameter :: linsat2 = 1.0, dturb_max = 100.0
58 integer,
parameter :: ener_norm =0
59 integer,
parameter :: ener_lsat=0
60 integer,
parameter :: nstdif = 1
61 integer,
parameter :: wave_sponge = 1
63 integer,
intent(in) :: levs
64 integer,
intent(in) :: im
65 integer,
intent(in) :: mpi_id, master, kdt
67 real(kind=kind_phys) ,
intent(in) :: dtp
68 real(kind=kind_phys) ,
intent(in) :: tau_ngw(im)
70 real(kind=kind_phys) ,
intent(in) :: vm(im,levs)
71 real(kind=kind_phys) ,
intent(in) :: um(im,levs)
72 real(kind=kind_phys) ,
intent(in) :: qm(im,levs)
73 real(kind=kind_phys) ,
intent(in) :: tm(im,levs)
75 real(kind=kind_phys) ,
intent(in) :: prsl(im,levs)
76 real(kind=kind_phys) ,
intent(in) :: prslk(im,levs)
77 real(kind=kind_phys) ,
intent(in) :: zmet(im,levs)
78 real(kind=kind_phys) ,
intent(in) :: prsi(im,levs+1)
79 real(kind=kind_phys) ,
intent(in) :: zmeti(im,levs+1)
80 real(kind=kind_phys) ,
intent(in) :: xlatd(im)
81 real(kind=kind_phys) ,
intent(in) :: sinlat(im)
82 real(kind=kind_phys) ,
intent(in) :: coslat(im)
86 real(kind=kind_phys) ,
intent(out) :: pdudt(im,levs)
87 real(kind=kind_phys) ,
intent(out) :: pdvdt(im,levs)
88 real(kind=kind_phys) ,
intent(out) :: pdtdt(im,levs)
89 real(kind=kind_phys) ,
intent(out) :: dked(im,levs)
90 real(kind=kind_phys) ,
intent(out) :: zngw(im)
96 real(kind=kind_phys) :: tauabs(im,levs)
97 real(kind=kind_phys) :: wrms(im,levs)
98 real(kind=kind_phys) :: trms(im,levs)
100 real(kind=kind_phys) :: zwrms(nwav,nazd), wrk1(levs), wrk2(levs)
101 real(kind=kind_phys) :: atrms(nazd, levs),awrms(nazd, levs), akzw(nwav,nazd, levs+1)
104 real(kind=kind_phys) :: taux(levs+1)
105 real(kind=kind_phys) :: tauy(levs+1)
106 real(kind=kind_phys) :: fpu(nazd, levs+1)
107 real(kind=kind_phys) :: ui(nazd, levs+1)
109 real(kind=kind_phys) :: fden_bn(levs+1)
110 real(kind=kind_phys) :: flux(nwav, nazd) , flux_m(nwav, nazd)
112 real(kind=kind_phys) :: bn(levs+1)
113 real(kind=kind_phys) :: bn2(levs+1)
114 real(kind=kind_phys) :: rhoint(levs+1)
115 real(kind=kind_phys) :: uint(levs+1)
116 real(kind=kind_phys) :: vint(levs+1)
117 real(kind=kind_phys) :: tint(levs+1)
119 real(kind=kind_phys) :: irhodz_mid(levs)
120 real(kind=kind_phys) :: suprf(levs+1)
121 real(kind=kind_phys) :: cstar(levs+1) ,cstar2(levs+1)
122 real(kind=kind_phys) :: v_zmet(levs+1)
123 real(kind=kind_phys) :: vueff(levs+1)
124 real(kind=kind_phys) :: dfdz_v(nazd, levs), dfdz_heat(nazd, levs)
126 real(kind=kind_phys),
dimension(levs) ::
atm , aum, avm, aqm, aprsl, azmet, dz_met
127 real(kind=kind_phys),
dimension(levs+1) :: aprsi, azmeti, dz_meti
129 real(kind=kind_phys),
dimension(levs) :: wrk3
130 real(kind=kind_phys),
dimension(levs) :: uold, vold, told, unew, vnew, tnew
131 real(kind=kind_phys),
dimension(levs) :: rho, rhomid, adif, cdif, acdif
132 real(kind=kind_phys),
dimension(levs) :: qmid, akt
133 real(kind=kind_phys),
dimension(levs+1) :: dktur, ktint, kvint
134 real(kind=kind_phys),
dimension(levs+1) :: fden_lsat, fden_bnen
136 integer,
dimension(levs) :: Anstab
138 real(kind=kind_phys) :: sig_u2az(nazd), sig_u2az_m(nazd)
139 real(kind=kind_phys) :: wave_dis(nwav, nazd), wave_disaz(nazd)
140 real(kind=kind_phys) :: rdci(nwav), rci(nwav)
141 real(kind=kind_phys) :: wave_act(nwav, nazd)
142 real(kind=kind_phys) :: ul(nazd)
146 real(kind=kind_phys) :: bvi, bvi2, bvi3, bvi4, rcms
147 real(kind=kind_phys) :: c2f2, cf1, wave_distot
150 real(kind=kind_phys) :: flux_norm
151 real(kind=kind_phys) :: taub_src, rho_src, zcool, vmdiff
153 real(kind=kind_phys) :: zthm, dtau, cgz, ucrit_maxdc
154 real(kind=kind_phys) :: vm_zflx_mode, vc_zflx_mode
155 real(kind=kind_phys) :: kzw2, kzw3, kdsat, cdf2, cdf1, wdop2,v_cdp2
156 real(kind=kind_phys) :: ucrit_max
157 real(kind=kind_phys) :: pwrms, ptrms
158 real(kind=kind_phys) :: zu, zcin, zcin2, zcin3, zcin4, zcinc
159 real(kind=kind_phys) :: zatmp, fluxs, zdep, ze1, ze2
162 real(kind=kind_phys) :: zdelp, zdelm, taud_min
163 real(kind=kind_phys) :: tvc, tvm, ptc, ptm
164 real(kind=kind_phys) :: umfp, umfm, umfc, ucrit3
165 real(kind=kind_phys) :: fmode, expdis, fdis
166 real(kind=kind_phys) :: v_kzi, v_kzw, v_cdp, v_wdp, tx1, fcorsat, dzcrit
167 real(kind=kind_phys) :: v_wdi, v_wdpc
168 real(kind=kind_phys) :: ugw, vgw, ek1, ek2, rdtp, rdtp2, rhp_wam
170 integer :: j, jj, k, kk, inc, jk, jkp, jl, iaz
171 integer :: ksrc, km2, km1, kp1, ktop
175 real(kind=kind_phys) :: uz, vz, shr2 , ritur, ktur
177 real(kind=kind_phys) :: kamp, zmetk, zgrow
178 real(kind=kind_phys) :: stab, stab_dt, dtstab
179 real(kind=kind_phys) :: nslope3
181 integer :: nstab, ist
182 real(kind=kind_phys) :: w1, w2, w3, dtdif
184 real(kind=kind_phys) :: dzmetm, dzmetp, dzmetf, bdif, bt_dif, apc, kturp
185 real(kind=kind_phys) :: rstar, rstar2
187 real(kind=kind_phys) :: snorm_ener, sigu2, flux_2_sig, ekin_norm
188 real(kind=kind_phys) :: taub_ch, sigu2_ch
189 real(kind=kind_phys) :: pr_kdis_eff, mf_diss_heat, ipr_max
190 real(kind=kind_phys) :: exp_sponge, mi_sponge, gipr
194 nslope3 = nslope + 3.0
195 pr_kdis_eff = gw_eff*pr_kdis
196 ipr_max = max(1.0, ipr_ktgw)
197 gipr = grav* ipr_ktgw
216 if (idebug_gwrms == 1)
then
217 tauabs=0.0; wrms =0.0 ; trms =0.0
226 ksrc= max(ilaunch, 3)
232 suprf(ktop) = kion(levs)
247 tx1 = omega2 * sinlat(j) *rv_kxw
250 ucrit_max = max(ucrit, cf1)
251 ucrit3 = ucrit_max*ucrit_max*ucrit_max
255 aprsl(1:levs) = prsl(jl,1:levs)
260 if (aprsl(k) .lt. psrc )
exit
262 ilaunch = max(k-1, 3)
263 ksrc= max(ilaunch, 3)
265 zngw(j) = zmet(j, ksrc)
273 aum(1:levs) = um(jl,1:levs)
274 avm(1:levs) = vm(jl,1:levs)
275 atm(1:levs) = tm(jl,1:levs)
276 aqm(1:levs) = qm(jl,1:levs)
277 azmet(1:levs) = zmet(jl,1:levs)
278 aprsi(1:levs+1) = prsi(jl,1:levs+1)
279 azmeti(1:levs+1) = zmeti(jl,1:levs+1)
281 rho_src = aprsl(ksrc)*rdi/
atm(ksrc)
282 taub_ch = max(tau_ngw(jl), tau_min)
286 sigu2 = taub_src/rho_src/v_kxw * zms
287 sig_u2az(1:nazd) = sigu2
292 dz_meti(jk) = azmeti(jk+1)-azmeti(jk)
293 dz_met(jk) = azmet(jk)-azmeti(jk-1)
299 tvc =
atm(jk)*(1. +fv*aqm(jk))
300 tvm =
atm(jk-1)*(1. +fv*aqm(jk-1))
301 ptc = tvc/ prslk(jl, jk)
302 ptm = tvm/prslk(jl,jk-1)
307 uint(jk) = 0.5*(aum(jk-1)+aum(jk))
308 vint(jk) = 0.5*(avm(jk-1)+avm(jk))
309 tint(jk) = 0.5*(tvc+tvm)
310 rhomid(jk) = aprsl(jk)*rdi/
atm(jk)
311 rhoint(jk) = aprsi(jk)*rdi*zthm
313 v_zmet(jk) = 2.*zdelp
314 zdelm = 1./dz_met(jk)
318 bn2(jk) = grav2cpd*zthm*(1.0+rcpdl*(tvc-tvm)*zdelm)
319 bn2(jk) = max(min(bn2(jk), bnv2max), bnv2min)
320 bn(jk) = sqrt(bn2(jk))
323 wrk3(jk)= 1./zdelp/rhomid(jk)
324 irhodz_mid(jk) = rdtp*zdelp*rhomid(jk)/rho_src
329 uz = aum(jk) - aum(jk-1)
330 vz = avm(jk) - avm(jk-1)
331 shr2 = (max(uz*uz+vz*vz, dw2min)) * zdelm *zdelm
333 zmetk = azmet(jk)* rh4
336 kamp = sqrt(shr2)*sc2 *zgrow
337 w1 = 1./(1. + 5*ritur)
338 ktur= min(max(kamp * w1 * w1, dked_min), dked_max)
339 zmetk = azmet(jk)* rhp
340 vueff(jk) = ktur + kvg(jk)
345 if (idebug_gwrms == 1)
then
347 wrk1(jk) = rv_kxw/rhoint(jk)
348 wrk2(jk)= rgrav2*zthm*zthm*bn2(jk)
357 rhoint(ktop) = 0.5*aprsi(levs)*rdi/
atm(jk)
358 tint(ktop) =
atm(jk)*(1. +fv*aqm(jk))
362 v_zmet(ktop) = v_zmet(jk)
363 vueff(ktop) = vueff(jk)
370 akt(jk) = -akt(jk)*(gor + (tint(jk+1)-tint(jk))/dz_meti(jk) )
374 bvi = bn(ksrc); bvi2 = bvi * bvi;
375 bvi3 = bvi2*bvi; bvi4 = bvi2 * bvi2; rcms = zms/bvi
380 ul(iaz) = zcosang(iaz) *uint(ksrc) + zsinang(iaz) *vint(ksrc)
385 cstar(jk) = bn(jk)/zms
386 cstar2(jk) = cstar(jk)*cstar(jk)
388 fden_lsat(jk) = rhoint(jk)/bn(jk)*v_kxw*linsat2
391 zu = zcosang(iaz)*uint(jk) + zsinang(iaz)*vint(jk)
396 rstar = 1./cstar(ksrc)
402 fpu(1:nazd, km2:ktop) =0.
406 zcin = zci(inc)*rstar
411 flux(inc,1) = rstar*(zcin*zcin)/(1.+ zcin**nslope3)
415 fpu(1,ksrc) = fpu(1,ksrc) + flux(inc,1)*zdci(inc)
418 akzw(inc, iaz, ksrc) = bvi*rci(inc)
425 flux_norm = taub_src / fpu(1, ksrc)
426 ze1 = flux_norm * bvi/rhoint(ksrc) *rstar *rstar2
428 fden_bn(jk) = ze1* rhoint(jk) / bn(jk)
432 flux(inc,1) = flux_norm*flux(inc,1)
436 if (ener_norm == 1)
then
439 zcin = zci(inc)*rstar
441 ze2 = zcin /(1.+ zcin**nslope3)
443 snorm_ener = snorm_ener + ze2*zdci(inc)*rstar
444 flux(inc,1) = ze2 * zcin
447 ekin_norm = 1./snorm_ener
453 ze1 = taub_src*zms/bvi * ekin_norm
457 flux(inc,1) = ze1* flux(inc,1)
458 taub_src = taub_src + flux(inc,1)*zdci(inc)
460 ze1 = ekin_norm * v_kxw * rstar2
462 fden_bnen(jk) = rhoint(jk) / bn(jk) *ze1
468 fpu(iaz, ksrc) = taub_src
469 fpu(iaz, km1) = taub_src
478 flux(inc,iaz) = flux(inc,1)
488 if (idebug_gwrms == 1)
then
491 tx1 = real(nazd)/rhoint(ksrc)*rv_kxw
495 ze1 = flux(inc,1)*zdci(inc)*tx1*v_kzw
497 ptrms = ptrms + ze1 * ze2
499 wrms(jl, ksrc) = pwrms
500 trms(jl, ksrc) = ptrms
512 sig_u2az_m(iaz) = sig_u2az(iaz)
516 umfc = .5*(umfm + umfp)
518 dfdz_v(iaz, jk) = 0.0
519 dfdz_heat(iaz, jk) = 0.0
525 flux_m(inc, iaz) = flux(inc, iaz)
530 if(wave_act(inc,iaz) == 1.0)
then
540 if (v_cdp .le. ucrit_max .or. cdf2 .le. 0.0)
then
544 wave_act(inc,iaz) =0.
545 akzw(inc, iaz, jkp) = pi/dz_meti(jk)
547 flux(inc,iaz) = fluxs
557 kzw2 = (bn2(jkp)-wdop2)/cdf2
561 if ( kzw2 > mkz2min )
then
563 akzw(inc, iaz, jkp) = v_kzw
571 v_wdp = v_kxw * v_cdp
572 v_wdi = kzw2*vueff(jk) + kion(jk)
573 v_wdpc = sqrt(v_wdp*v_wdp +v_wdi*v_wdi)
574 v_kzi = v_kzw*v_wdi/v_wdpc
577 ze1 = v_kzi*v_zmet(jk)
579 if (ze1 .ge. 1.e-2)
then
580 expdis = max(exp(-ze1), 0.01)
582 expdis = 1./(1.+ ze1)
586 wave_act(inc,iaz) = 1.0
587 fmode = flux(inc,iaz)
589 flux_2_sig = v_kzw/v_kxw/rhoint(jkp)
590 w1 = v_wdpc/kzw2/v_kzw/v_zmet(jk)
597 wave_act(inc,iaz) = 0.0
598 akzw(inc, iaz, jkp) = v_kzw
604 fdis = fmode*expdis*wave_act(inc,iaz)
627 if (ener_norm == 0) fluxs= fden_bn(jkp)*cdf2*wave_act(inc,iaz)
631 if (ener_lsat == 1) fluxs= fden_lsat(jkp)*cdf2*sqrt(cdf2)*rdci(inc)*wave_act(inc,iaz)
633 if (ener_norm == 1)
then
637 if (ener_lsat == 0) fluxs= fden_bnen(jk)*cdf2*wave_act(inc,iaz)*sig_u2az_m(iaz)
641 if (ener_lsat == 1) fluxs= fden_lsat(jkp)*cdf2*sqrt(cdf2)*rdci(inc)*wave_act(inc,iaz)
653 flux(inc,iaz) = fluxs
654 ze2 = log(ze1/fluxs)*w1
661 dtau = flux_m(inc,iaz)-flux(inc,iaz)
662 if (dtau .lt. 0)
then
663 flux(inc,iaz) = flux_m(inc,iaz)
668 if ( azmeti(jkp) .ge. zsp_gw)
then
669 mi_sponge = .5/dz_meti(jk)
670 ze2 = v_wdp /v_kzw * mi_sponge
671 v_wdi = ze2 + v_wdi*0.25
672 v_wdpc = sqrt(v_wdp*v_wdp +v_wdi*v_wdi)
673 v_kzi = v_kzw*v_wdi/v_wdpc
675 ze1 = v_kzi*v_zmet(jk)
676 exp_sponge = exp(-ze1)
680 flux(inc,iaz) = flux(inc,iaz) *exp_sponge
688 if (wave_act(inc,iaz) == 1)
then
691 vc_zflx_mode = flux(inc,iaz)
692 vmdiff = max(0., flux_m(inc,iaz)-vc_zflx_mode)
693 if (vmdiff <= 0. ) vc_zflx_mode = flux_m(inc,iaz)
694 ze1 = vc_zflx_mode*zcinc
695 fpu(iaz, jkp) = fpu(iaz,jkp) + ze1
696 sig_u2az(iaz) = sig_u2az(iaz) + ze1*flux_2_sig
705 zdelp = wrk3(jk)* v_cdp *zcinc * vmdiff
714 dfdz_v(iaz, jk) = dfdz_v(iaz,jk) + zdelp
715 dfdz_heat(iaz, jk) = dfdz_heat(iaz,jk) + zdelp
721 if (fpu(iaz, jkp) > ze1 ) fpu(iaz, jkp) = ze1
725 if (idebug_gwrms == 1)
then
729 if (wave_act(inc,iaz) > 0.)
then
730 v_kzw =akzw(inc, iaz, jk)
731 ze1 = flux(inc,iaz)*v_kzw*zdci(inc)*wrk1(jk)
733 ptrms = ptrms + ze1*wrk2(jk)
736 awrms(iaz, jk) = pwrms
737 atrms(iaz, jk) = ptrms
746 tx1 = sum(abs(dfdz_heat(1:nazd, jk)))/bn2(jk)
747 ze1=max(dked_min, tx1)
748 ze2=min(dked_max, ze1)
749 vueff(jkp) = ze2 + vueff(jkp)
771 taux(jk) = taux(jk) + fpu(iaz,jk)*zcosang(iaz)
772 tauy(jk) = tauy(jk) + fpu(iaz,jk)*zsinang(iaz)
773 pdtdt(jl,jk) = pdtdt(jl,jk) + dfdz_v(iaz,jk)
774 dked(jl,jk) = dked(jl,jk) + dfdz_heat(iaz,jk)
777 jk = ktop; taux(jk)=0.; tauy(jk)=0.
779 taux(jk) = taux(jk) + fpu(iaz,jk)*zcosang(iaz)
780 tauy(jk) = tauy(jk) + fpu(iaz,jk)*zsinang(iaz)
783 if (idebug_gwrms == 1)
then
786 wrms(jl,jk) =wrms(jl,jk) + awrms(iaz,jk)
787 trms(jl,jk) =trms(jl,jk) + atrms(iaz,jk)
788 tauabs(jl,jk)=tauabs(jl,jk) + fpu(iaz,jk)
796 zdelp = wrk3(jk)*gw_eff
797 ze1 = (taux(jkp)-taux(jk))* zdelp
798 ze2 = (tauy(jkp)-tauy(jk))* zdelp
800 if (abs(ze1) >= maxdudt )
then
801 ze1 = sign(maxdudt, ze1)
803 if (abs(ze2) >= maxdudt )
then
804 ze2 = sign(maxdudt, ze2)
813 if (knob_ugwp_doheat == 1)
then
817 pdtdt(jl,jk) = pdtdt(jl,jk)*gw_eff
819 if (abs(ze2) >= max_eps ) pdtdt(jl,jk) = sign(max_eps, ze2)
821 dked(jl,jk) = dked(jl,jk)/bn2(jk)
822 ze1 = max(dked_min, dked(jl,jk))
823 dked(jl,jk) = min(dked_max, ze1)
824 qmid(jk) = pdtdt(j,jk)
834 dktur(1:levs) = dked(jl,1:levs)
838 adif(jk) =.25*(dktur(jk-1)+ dktur(jk+1)) + .5*dktur(jk)
840 dktur(ksrc:levs-1) = adif(ksrc:levs-1)
842 dktur(levs) = .5*( dked(jl,levs)+ dked(jl,levs-1))
843 dktur(levs+1) = dktur(levs)
846 ze1 = .5*( dktur(jk) +dktur(jk-1) )
848 ktint(jk) = ze1*ipr_ktgw
855 ze2 = qmid(jk) + dktur(jk)*akt(jk) + grav*(ktint(jk+1)-ktint(jk))/dz_meti(jk)
857 if (abs(ze2) >= max_eps ) qmid(jk) = sign(max_eps, ze2)
858 pdtdt(jl,jk) = qmid(jk)*rcpd
859 dked(jl, jk) = dktur(jk)
868 uold(km2:levs) = aum(km2:levs)+pdudt(jl,km2:levs)*dtp
869 vold(km2:levs) = avm(km2:levs)+pdvdt(jl,km2:levs)*dtp
870 told(km2:levs) =
atm(km2:levs)+pdtdt(jl,km2:levs)*dtp
875 dktur(km2:levs) = dked_min
878 uz = uold(jk) - uold(jk-1)
879 vz = vold(jk) - vold(jk-1)
883 tvc = told(jk) * (1. +fv*aqm(jk))
884 tvm = told(jk-1) * (1. +fv*aqm(jk-1))
885 zthm = 2.0 / (tvc+tvm)
886 shr2 = (max(uz*uz+vz*vz, dw2min)) * zdelm *zdelm
888 bn2(jk) = grav2cpd*zthm * (1.0+rcpdl*(tvc-tvm)*zdelm)
890 bn2(jk) = max(min(bn2(jk), bnv2max), bnv2min)
891 zmetk = azmet(jk)* rh4
894 w1 = 1./(1. + 5*ritur)
895 ze2 = min( sc2 *zgrow, 4.*ze1*ze1)
899 kamp = sqrt(shr2)* ze2 * w1 * w1
900 ktur= min(max(kamp, dked_min), dked_max)
905 dked(jl, jk) = dked(jl, jk) +ktur
912 if (knob_ugwp_dokdis == 2)
then
915 ze1 = min(.5*(dktur(jk) +dktur(jk-1)), dturb_max)
916 kvint(jk) = kvint(jk) + ze1
919 kvint(km1) = kvint(ksrc)
920 kvint(ktop) = kvint(levs)
922 dzmetm = 1./dz_met(km1)
927 dzmetp = 1./dz_met(jk+1)
928 dzmetf = 1./(dz_meti(jk)*rhomid(jk))
931 ktur = kvint(jk) *rhoint(jk) * dzmetf
932 kturp =kvint(jk+1)*rhoint(jk+1) * dzmetf
934 adif(jk) = ktur * dzmetm
935 cdif(jk) = kturp * dzmetp
936 apc = adif(jk)+cdif(jk)
941 anstab(jk) = floor(w1*dtp) + 1
948 nstab = maxval( anstab(ksrc:levs-1))
954 dtdif = dtp/real(nstab)
959 bdif = ze1 - acdif(k)
960 bt_dif = ze1 - acdif(k)* ipr_ktgw
961 unew(k) = uold(k)*bdif + uold(k-1)*adif(k) + uold(k+1)*cdif(k)
962 vnew(k) = vold(k)*bdif + vold(k-1)*adif(k) + vold(k+1)*cdif(k)
963 tnew(k) = told(k)*bt_dif+(told(k-1)*adif(k) + told(k+1)*cdif(k))*ipr_ktgw
966 uold(ksrc:levs-1) = unew(ksrc:levs-1)*dtdif
967 vold(ksrc:levs-1) = vnew(ksrc:levs-1)*dtdif
968 told(ksrc:levs-1) = tnew(ksrc:levs-1)*dtdif
972 uold(levs) = uold(levs-1)
973 vold(levs) = vold(levs-1)
974 told(levs) = told(levs-1)
983 ze2 = rdtp*(uold(k) - aum(k))
984 ze1 = rdtp*(vold(k) - avm(k))
985 pdtdt(jl,k)= rdtp*( told(k) -
atm(k) )
987 if (abs(pdtdt(jl,k)) >= maxdtdt ) pdtdt(jl,k) = sign(maxdtdt,pdtdt(jl,k) )
988 if (abs(ze1) >= maxdudt )
then
989 ze1 = sign(maxdudt, ze1)
991 if (abs(ze2) >= maxdudt )
then
992 ze2 = sign(maxdudt, ze2)
997 uz = uold(k+1) - uold(k-1)
998 vz = vold(k+1) - vold(k-1)
999 ze2 = 1./(dz_met(k+1)+dz_met(k) )
1000 mf_diss_heat = rcpd*kvint(k)*(uz*uz +vz*vz)*ze2*ze2
1001 pdtdt(jl,k)= pdtdt(jl,k) + mf_diss_heat
1013 if (kdt ==1 .and. mpi_id == master)
then
1015 print *,
' ugwpv1: nazd-nw-ilaunch=', nazd, nwav,ilaunch, maxval(kvg),
' kvg '
1016 print *,
'ugwpv1: zdci(inc)=' , maxval(zdci), minval(zdci)
1017 print *,
'ugwpv1: zcimax=' , maxval(zci) ,
' zcimin=' , minval(zci)
1024 if (kdt == 1 .and. mpi_id == master)
then
1025 print *,
'vgw done nstab ', nstab
1027 print *, maxval(pdudt)*86400., minval(pdudt)*86400,
'vgw ax ugwp'
1028 print *, maxval(pdvdt)*86400., minval(pdvdt)*86400,
'vgw ay ugwp'
1029 print *, maxval(dked)*1., minval(dked)*1,
'vgw keddy m2/sec ugwp'
1030 print *, maxval(pdtdt)*86400., minval(pdtdt)*86400,
'vgw eps ugwp'