subroutine turbdiff(ngrid,nlay,nq,rnat, & ptimestep,pcapcal,lecrit, & pplay,pplev,pzlay,pzlev,pz0, & pu,pv,pt,ppopsk,pq,ptsrf,pemis,pqsurf, & pdufi,pdvfi,pdtfi,pdqfi,pfluxsrf, & Pdudif,pdvdif,pdtdif,pdtsrf,sensibFlux,pq2, & pdqdif,pdqsdif,lastcall) use watercommon_h, only : RLVTT, To, RCPD, mx_eau_sol use radcommon_h, only : sigma implicit none !================================================================== ! ! Purpose ! ------- ! Turbulent diffusion (mixing) for pot. T, U, V and tracers ! ! Implicit scheme ! We start by adding to variables x the physical tendencies ! already computed. We resolve the equation: ! ! x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1) ! ! Authors ! ------- ! F. Hourdin, F. Forget, R. Fournier (199X) ! R. Wordsworth, B. Charnay (2010) ! J. Leconte (2012): To f90 ! - Rewritten the diffusion scheme to conserve total enthalpy ! by accounting for dissipation of turbulent kinetic energy. ! - Accounting for lost mean flow kinetic energy should come soon. ! !================================================================== !----------------------------------------------------------------------- ! declarations ! ------------ #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "callkeys.h" #include "surfdat.h" #include "comgeomfi.h" #include "tracer.h" #include "watercap.h" ! arguments ! --------- INTEGER ngrid,nlay REAL ptimestep REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1) REAL pu(ngrid,nlay),pv(ngrid,nlay) REAL pt(ngrid,nlay),ppopsk(ngrid,nlay) REAL ptsrf(ngrid),pemis(ngrid) REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay) REAL pdtfi(ngrid,nlay) REAL pfluxsrf(ngrid) REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay) REAL pdtdif(ngrid,nlay) REAL pdtsrf(ngrid),sensibFlux(ngrid),pcapcal(ngrid) REAL pq2(ngrid,nlay+1) integer rnat(ngrid) ! Arguments added for condensation logical lecrit REAL pz0 ! Tracers ! -------- integer nq real pqsurf(ngrid,nq) real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) real pdqdif(ngrid,nlay,nq) real pdqsdif(ngrid,nq) ! local ! ----- integer ilev,ig,ilay,nlev REAL z4st,zdplanck(ngridmx) REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1) REAL zcdv(ngridmx),zcdh(ngridmx) REAL zcdv_true(ngridmx),zcdh_true(ngridmx) REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx) REAL zh(ngridmx,nlayermx),zt(ngridmx,nlayermx) REAL ztsrf(ngridmx) REAL z1(ngridmx),z2(ngridmx) REAL zmass(ngridmx,nlayermx) REAL zfluxv(ngridmx,nlayermx),zfluxt(ngridmx,nlayermx),zfluxq(ngridmx,nlayermx) REAL zb0(ngridmx,nlayermx) REAL zExner(ngridmx,nlayermx),zovExner(ngridmx,nlayermx) REAL zcv(ngridmx,nlayermx),zdv(ngridmx,nlayermx) !inversion coefficient for winds REAL zct(ngridmx,nlayermx),zdt(ngridmx,nlayermx) !inversion coefficient for temperature REAL zcq(ngridmx,nlayermx),zdq(ngridmx,nlayermx) !inversion coefficient for tracers REAL zcst1 REAL zu2!, a REAL zcq0(ngridmx),zdq0(ngridmx) REAL zx_alf1(ngridmx),zx_alf2(ngridmx) LOGICAL firstcall SAVE firstcall LOGICAL lastcall ! Tracers ! ------- INTEGER iq REAL zq(ngridmx,nlayermx,nqmx) REAL rho(ngridmx) ! near-surface air density REAL qsat(ngridmx) DATA firstcall/.true./ REAL kmixmin ! Variables added for implicit latent heat inclusion ! -------------------------------------------------- real dqsat(ngridmx), qsat_temp1, qsat_temp2 integer ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface... save ivap, iliq, iliq_surf,iice_surf real, parameter :: karman=0.4 real cd0, roughratio real dqsdif_total(ngrid) real zq0(ngrid) ! Coherence test ! -------------- IF (firstcall) THEN if(water)then ivap=igcm_h2o_vap iliq=igcm_h2o_ice iliq_surf=igcm_h2o_vap iice_surf=igcm_h2o_ice ! simply to make the code legible ! to be generalised later endif sensibFlux(:)=0. firstcall=.false. ENDIF !----------------------------------------------------------------------- ! 1. Initialisation ! ----------------- nlev=nlay+1 ! Calculate rho*dz, (P/Ps)**(R/cp) and dt*rho/dz=dt*rho**2 g/dp ! with rho=p/RT=p/ (R Theta) (p/ps)**kappa ! --------------------------------------------- DO ilay=1,nlay DO ig=1,ngrid zmass(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g zExner(ig,ilay)=(pplev(ig,ilay)/pplev(ig,1))**rcp zovExner(ig,ilay)=1./ppopsk(ig,ilay) ENDDO ENDDO zcst1=4.*g*ptimestep/(R*R) DO ilev=2,nlev-1 DO ig=1,ngrid zb0(ig,ilev)=pplev(ig,ilev)/(pt(ig,ilev-1)+pt(ig,ilev)) zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/(pplay(ig,ilev-1)-pplay(ig,ilev)) ENDDO ENDDO DO ig=1,ngrid zb0(ig,1)=ptimestep*pplev(ig,1)/(R*ptsrf(ig)) ENDDO dqsdif_total(:)=0.0 !----------------------------------------------------------------------- ! 2. Add the physical tendencies computed so far ! ---------------------------------------------- DO ilev=1,nlay DO ig=1,ngrid zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep zh(ig,ilev)=pt(ig,ilev)*zovExner(ig,ilev) !for call vdif_kc, but could be moved and computed there ENDDO ENDDO if(tracer) then DO iq =1, nq DO ilev=1,nlay DO ig=1,ngrid zq(ig,ilev,iq)=pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep ENDDO ENDDO ENDDO end if !----------------------------------------------------------------------- ! 3. Turbulence scheme ! -------------------- ! ! Source of turbulent kinetic energy at the surface ! ------------------------------------------------- ! Formula is Cd_0 = (karman / log[1+z1/z0])^2 DO ig=1,ngrid roughratio = 1. + pzlay(ig,1)/pz0 cd0 = karman/log(roughratio) cd0 = cd0*cd0 zcdv_true(ig) = cd0 zcdh_true(ig) = cd0 ! zcdv_true(ig)=0.!JL12 uncomment to disable atm/surface momentum flux ! zcdh_true(ig)=0.!JL12 uncomment to disable sensible heat flux ENDDO DO ig=1,ngrid zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) zcdv(ig)=zcdv_true(ig)*sqrt(zu2) zcdh(ig)=zcdh_true(ig)*sqrt(zu2) ENDDO ! Turbulent diffusion coefficients in the boundary layer ! ------------------------------------------------------ call vdif_kc(ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature ! Adding eddy mixing to mimic 3D general circulation in 1D ! R. Wordsworth & F. Forget (2010) if ((ngrid.eq.1)) then kmixmin = 1.0e-2 ! minimum eddy mix coeff in 1D do ilev=1,nlay do ig=1,ngrid zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev)) zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev)) end do end do end if !JL12 change zkv at the surface by zcdv to calculate the surface momentum properly DO ig=1,ngrid zkv(ig,1)=zcdv(ig) ENDDO !we treat only winds, energy and tracers coefficients will be computed with upadted winds !JL12 calculate the flux coefficients (tables multiplied elements by elements) zfluxv=zkv(:,1:nlay)*zb0 !----------------------------------------------------------------------- ! 4. Implicit inversion of u ! -------------------------- ! u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1) ! avec ! /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.) ! et ! dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1) ! donc les entrees sont /zcdv/ pour la condition a la limite sol ! et /zkv/ = Ku DO ig=1,ngrid z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay)) zcv(ig,nlay)=zmass(ig,nlay)*zu(ig,nlay)*z1(ig) zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig) ENDDO DO ilay=nlay-1,1,-1 DO ig=1,ngrid z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay) + zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1))) zcv(ig,ilay)=(zmass(ig,ilay)*zu(ig,ilay)+zfluxv(ig,ilay+1)*zcv(ig,ilay+1))*z1(ig) zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig) ENDDO ENDDO DO ig=1,ngrid zu(ig,1)=zcv(ig,1) ENDDO DO ilay=2,nlay DO ig=1,ngrid zu(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zu(ig,ilay-1) ENDDO ENDDO !----------------------------------------------------------------------- ! 5. Implicit inversion of v ! -------------------------- ! v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1) ! avec ! /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.) ! et ! dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1) ! donc les entrees sont /zcdv/ pour la condition a la limite sol ! et /zkv/ = Kv DO ig=1,ngrid z1(ig)=1./(zmass(ig,nlay)+zfluxv(ig,nlay)) zcv(ig,nlay)=zmass(ig,nlay)*zv(ig,nlay)*z1(ig) zdv(ig,nlay)=zfluxv(ig,nlay)*z1(ig) ENDDO DO ilay=nlay-1,1,-1 DO ig=1,ngrid z1(ig)=1./(zmass(ig,ilay)+zfluxv(ig,ilay)+zfluxv(ig,ilay+1)*(1.-zdv(ig,ilay+1))) zcv(ig,ilay)=(zmass(ig,ilay)*zv(ig,ilay)+zfluxv(ig,ilay+1)*zcv(ig,ilay+1))*z1(ig) zdv(ig,ilay)=zfluxv(ig,ilay)*z1(ig) ENDDO ENDDO DO ig=1,ngrid zv(ig,1)=zcv(ig,1) ENDDO DO ilay=2,nlay DO ig=1,ngrid zv(ig,ilay)=zcv(ig,ilay)+zdv(ig,ilay)*zv(ig,ilay-1) ENDDO ENDDO !---------------------------------------------------------------------------- ! 6. Implicit inversion of h, not forgetting the coupling with the ground ! h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1) ! avec ! /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.) ! et ! dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1) ! donc les entrees sont /zcdh/ pour la condition de raccord au sol ! et /zkh/ = Kh ! Using the wind modified by friction for lifting and sublimation ! --------------------------------------------------------------- DO ig=1,ngrid zu2 = zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) zcdv(ig) = zcdv_true(ig)*sqrt(zu2) zcdh(ig) = zcdh_true(ig)*sqrt(zu2) zkh(ig,1)=zcdh(ig) ENDDO ! JL12 calculate the flux coefficients (tables multiplied elements by elements) ! --------------------------------------------------------------- zfluxq=zkh(:,1:nlay)*zb0 !JL12 we save zfluxq which doesn't need the Exner factor zfluxt=zfluxq*zExner DO ig=1,ngrid z1(ig)=1./(zmass(ig,nlay)+zfluxt(ig,nlay)*zovExner(ig,nlay)) zct(ig,nlay)=zmass(ig,nlay)*zt(ig,nlay)*z1(ig) zdt(ig,nlay)=zfluxt(ig,nlay)*zovExner(ig,nlay-1)*z1(ig) ENDDO DO ilay=nlay-1,2,-1 DO ig=1,ngrid z1(ig)=1./(zmass(ig,ilay)+zfluxt(ig,ilay)*zovExner(ig,ilay)+ & zfluxt(ig,ilay+1)*(zovExner(ig,ilay)-zdt(ig,ilay+1)*zovExner(ig,ilay+1))) zct(ig,ilay)=(zmass(ig,ilay)*zt(ig,ilay)+zfluxt(ig,ilay+1)*zct(ig,ilay+1)*zovExner(ig,ilay+1))*z1(ig) zdt(ig,ilay)=zfluxt(ig,ilay)*z1(ig)*zovExner(ig,ilay-1) ENDDO ENDDO !JL12 we treat last point afterward because zovExner(ig,ilay-1) does not exist there DO ig=1,ngrid z1(ig)=1./(zmass(ig,1)+zfluxt(ig,1)*zovExner(ig,1)+ & zfluxt(ig,2)*(zovExner(ig,1)-zdt(ig,2)*zovExner(ig,2))) zct(ig,1)=(zmass(ig,1)*zt(ig,1)+zfluxt(ig,2)*zct(ig,2)*zovExner(ig,2))*z1(ig) zdt(ig,1)=zfluxt(ig,1)*z1(ig) ENDDO ! Calculate (d Planck / dT) at the interface temperature ! ------------------------------------------------------ z4st=4.0*sigma*ptimestep DO ig=1,ngrid zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) ENDDO ! Calculate temperature tendency at the interface (dry case) ! ---------------------------------------------------------- ! Sum of fluxes at interface at time t + \delta t gives change in T: ! radiative fluxes ! turbulent convective (sensible) heat flux ! flux (if any) from subsurface if(.not.water) then DO ig=1,ngrid z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) & + pfluxsrf(ig)*ptimestep + zdplanck(ig)*ptsrf(ig) z2(ig) = pcapcal(ig)+zdplanck(ig)+cpp*zfluxt(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) ztsrf(ig) = z1(ig) / z2(ig) pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig) ENDDO ! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme ! Recalculate temperature to top of atmosphere, starting from ground ! ------------------------------------------------------------------ DO ilay=2,nlay DO ig=1,ngrid zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1) ENDDO ENDDO endif ! not water !----------------------------------------------------------------------- ! TRACERS (no vapour) ! ------- if(tracer) then ! Calculate vertical flux from the bottom to the first layer (dust) ! ----------------------------------------------------------------- do ig=1,ngridmx rho(ig) = zb0(ig,1) /ptimestep end do pdqsdif(:,:)=0. ! Implicit inversion of q ! ----------------------- do iq=1,nq if (iq.ne.ivap) then DO ig=1,ngrid z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay)) zcq(ig,nlay)=zmass(ig,nlay)*zq(ig,nlay,iq)*z1(ig) zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig) ENDDO DO ilay=nlay-1,2,-1 DO ig=1,ngrid z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1))) zcq(ig,ilay)=(zmass(ig,ilay)*zq(ig,ilay,iq)+zfluxq(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig) zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig) ENDDO ENDDO if ((water).and.(iq.eq.iliq)) then ! special case for condensed water tracer: do not include ! h2o ice tracer from surface (which is set when handling ! h2o vapour case (see further down). ! zb(ig,1)=0 if iq ne ivap DO ig=1,ngrid z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2))) zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2))*z1(ig) ENDDO else ! general case DO ig=1,ngrid z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2))) zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2)+(-pdqsdif(ig,iq))*ptimestep)*z1(ig) ! tracer flux from surface ! currently pdqsdif always zero here, ! so last line is superfluous enddo endif ! of if (water.and.(iq.eq.igcm_h2o_ice)) ! Starting upward calculations for simple tracer mixing (e.g., dust) do ig=1,ngrid zq(ig,1,iq)=zcq(ig,1) end do do ilay=2,nlay do ig=1,ngrid zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq) end do end do endif ! if (iq.ne.ivap) ! Calculate temperature tendency including latent heat term ! and assuming an infinite source of water on the ground ! ------------------------------------------------------------------ if (water.and.(iq.eq.ivap)) then ! compute evaporation efficiency do ig = 1, ngrid if(rnat(ig).eq.1)then dryness(ig)=pqsurf(ig,iliq_surf)+pqsurf(ig,iice_surf) dryness(ig)=MIN(1.,2*dryness(ig)/mx_eau_sol) dryness(ig)=MAX(0.,dryness(ig)) endif enddo do ig=1,ngrid ! Calculate the value of qsat at the surface (water) call watersat(ptsrf(ig),pplev(ig,1),qsat(ig)) call watersat(ptsrf(ig)-0.0001,pplev(ig,1),qsat_temp1) call watersat(ptsrf(ig)+0.0001,pplev(ig,1),qsat_temp2) dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002 ! calculate dQsat / dT by finite differences ! we cannot use the updated temperature value yet... enddo ! coefficients for q do ig=1,ngrid z1(ig)=1./(zmass(ig,nlay)+zfluxq(ig,nlay)) zcq(ig,nlay)=zmass(ig,nlay)*zq(ig,nlay,iq)*z1(ig) zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig) enddo do ilay=nlay-1,2,-1 do ig=1,ngrid z1(ig)=1./(zmass(ig,ilay)+zfluxq(ig,ilay)+zfluxq(ig,ilay+1)*(1.-zdq(ig,ilay+1))) zcq(ig,ilay)=(zmass(ig,ilay)*zq(ig,ilay,iq)+zfluxq(ig,ilay+1)*zcq(ig,ilay+1))*z1(ig) zdq(ig,ilay)=zfluxq(ig,ilay)*z1(ig) enddo enddo do ig=1,ngrid z1(ig)=1./(zmass(ig,1)+zfluxq(ig,1)*dryness(ig)+zfluxq(ig,2)*(1.-zdq(ig,2))) zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2))*z1(ig) zdq(ig,1)=dryness(ig)*zfluxq(ig,1)*z1(ig) enddo do ig=1,ngrid !calculation of surface temperature zdq0(ig) = dqsat(ig) zcq0(ig) = qsat(ig)-dqsat(ig)*ptsrf(ig) z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1) & + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep & + zfluxq(ig,1)*dryness(ig)*RLVTT*((zdq(ig,1)-1.0)*zcq0(ig)+zcq(ig,1)) z2(ig) = pcapcal(ig) + cpp*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) & + zdplanck(ig)+zfluxq(ig,1)*dryness(ig)*RLVTT*zdq0(ig)*(1.0-zdq(ig,1)) ztsrf(ig) = z1(ig) / z2(ig) ! calculation of qs and q1 zq0(ig) = zcq0(ig)+zdq0(ig)*ztsrf(ig) zq(ig,1,iq) = zcq(ig,1)+zdq(ig,1)*zq0(ig) ! calculation of evaporation dqsdif_total(ig)=zfluxq(ig,1)*dryness(ig)*(zq(ig,1,ivap)-zq0(ig)) ! -------------------------------------------------------- ! Now check if we've taken too much water from the surface ! This can only occur on the continent ! If we do, we recompute Tsurf, T1 and q1 accordingly if((-dqsdif_total(ig).gt.(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf))).and.rnat(ig).eq.1)then !water flux * ptimestep dqsdif_total(ig)=-(pqsurf(ig,iice_surf)+pqsurf(ig,iliq_surf)) !recompute surface temperature z1(ig) = pcapcal(ig)*ptsrf(ig) +cpp*zfluxq(ig,1)*zct(ig,1)*zovExner(ig,1) & + zdplanck(ig)*ptsrf(ig) + pfluxsrf(ig)*ptimestep & + RLVTT*dqsdif_total(ig) z2(ig) = pcapcal(ig) + cpp*zfluxq(ig,1)*(1.-zovExner(ig,1)*zdt(ig,1)) & + zdplanck(ig) ztsrf(ig) = z1(ig) / z2(ig) !recompute q1 with new water flux from surface zq(ig,1,iq) = (zmass(ig,1)*(pq(ig,1,iq)+ptimestep*pdqfi(ig,1,iq)) & +zfluxq(ig,2)*zcq(ig,1)-dqsdif_total(ig)) & / (zmass(ig,1)+(1.-zdq(ig,2))*zfluxq(ig,2)) end if ! calculation surface T tendency and T(1) pdtsrf(ig) = (ztsrf(ig) - ptsrf(ig))/ptimestep zt(ig,1) = zct(ig,1) + zdt(ig,1)*ztsrf(ig) enddo ! recalculate temperature and q(vap) to top of atmosphere, starting from ground do ilay=2,nlay do ig=1,ngrid zq(ig,ilay,iq)=zcq(ig,ilay)+zdq(ig,ilay)*zq(ig,ilay-1,iq) zt(ig,ilay)=zct(ig,ilay)+zdt(ig,ilay)*zt(ig,ilay-1) end do end do do ig=1,ngrid ! -------------------------------------------------------------------------- ! On the ocean, if T > 0 C then the vapour tendency must replace the ice one ! The surface vapour tracer is actually liquid. To make things difficult. if (rnat(ig).eq.0) then ! unfrozen ocean pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep pdqsdif(ig,iice_surf)=0.0 elseif (rnat(ig).eq.1) then ! (continent) ! If water is evaporating / subliming, we take it from ice before liquid ! -- is this valid?? if(dqsdif_total(ig).lt.0)then if (-dqsdif_total(ig).gt.pqsurf(ig,iice_surf))then pdqsdif(ig,iice_surf) = -pqsurf(ig,iice_surf)/ptimestep ! removes all the ice! pdqsdif(ig,iliq_surf) = dqsdif_total(ig)/ptimestep- pdqsdif(ig,iice_surf) ! take the remainder from the liquid instead else pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep pdqsdif(ig,iliq_surf)=0. end if else !dqsdif_total(ig).ge.0 !If water vapour is condensing, we must decide whether it forms ice or liquid. if(ztsrf(ig).gt.To)then pdqsdif(ig,iice_surf)=0.0 pdqsdif(ig,iliq_surf)=dqsdif_total(ig)/ptimestep else pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep pdqsdif(ig,iliq_surf)=0.0 endif endif elseif (rnat(ig).eq.2) then ! (continental glaciers) pdqsdif(ig,iliq_surf)=0.0 pdqsdif(ig,iice_surf)=dqsdif_total(ig)/ptimestep endif !rnat end do ! of DO ig=1,ngrid endif ! if (water et iq=ivap) end do ! of do iq=1,nq endif ! traceur !----------------------------------------------------------------------- ! 8. Final calculation of the vertical diffusion tendencies ! ----------------------------------------------------------------- do ilev = 1, nlay do ig=1,ngrid pdudif(ig,ilev)=(zu(ig,ilev)-(pu(ig,ilev)+pdufi(ig,ilev)*ptimestep))/ptimestep pdvdif(ig,ilev)=(zv(ig,ilev)-(pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep))/ptimestep pdtdif(ig,ilev)=( zt(ig,ilev)- pt(ig,ilev))/ptimestep-pdtfi(ig,ilev) enddo enddo DO ig=1,ngrid ! computing sensible heat flux (atm => surface) sensibFlux(ig)=cpp*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig)) ENDDO if (tracer) then do iq = 1, nq do ilev = 1, nlay do ig=1,ngrid pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-(pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep))/ptimestep enddo enddo enddo endif if(water)then call writediagfi(ngrid,'beta','Dryness coefficient',' ',2,dryness) endif ! if(lastcall)then ! if(ngrid.eq.1)then ! print*,'Saving k.out...' ! OPEN(12,file='k.out',form='formatted') ! DO ilay=1,nlay ! write(12,*) zkh(1,ilay), pplay(1,ilay) ! ENDDO ! CLOSE(12) ! endif ! endif return end