Index: DZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F (revision 1471) +++ (revision ) @@ -1,194 +1,0 @@ -! -! $Id$ -! - SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) - -c Auteur : P. Le Van . -c - IMPLICIT NONE - -#include "dimensions.h" -#include "paramet.h" -#include "iniprint.h" -#include "logic.h" -c -c======================================================================= -c -c -c s = sigma ** kappa : coordonnee verticale -c dsig(l) : epaisseur de la couche l ds la coord. s -c sig(l) : sigma a l'interface des couches l et l-1 -c ds(l) : distance entre les couches l et l-1 en coord.s -c -c======================================================================= -c - REAL pa,preff - REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1) - REAL presnivs(llm) -c -c declarations: -c ------------- -c - REAL sig(llm+1),dsig(llm) - real zzz(1:llm+1) - real dzz(1:llm) - real zk,zkm1,dzk1,dzk2,k0,k1 -c - INTEGER l - REAL snorm,dsigmin - REAL alpha,beta,gama,delta,deltaz,h - INTEGER np,ierr - REAL pi,x - - REAL SSUM -c -c----------------------------------------------------------------------- -c - pi=2.*ASIN(1.) - - OPEN(99,file='sigma.def',status='old',form='formatted', - s iostat=ierr) - -c----------------------------------------------------------------------- -c cas 1 on lit les options dans sigma.def: -c ---------------------------------------- - - IF (ierr.eq.0) THEN - - READ(99,*) h ! hauteur d'echelle 8. - READ(99,*) deltaz ! epaiseur de la premiere couche 0.04 - READ(99,*) beta ! facteur d'acroissement en haut 1.3 - READ(99,*) k0 ! nombre de couches dans la transition surf - READ(99,*) k1 ! nombre de couches dans la transition haute - CLOSE(99) - alpha=deltaz/(llm*h) - write(lunout,*)'h,alpha,k0,k1,beta' - -c read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2 - - alpha=deltaz/tanh(1./k0)*2. - zkm1=0. - sig(1)=1. - do l=1,llm - sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) - + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) - zk=-h*log(sig(l+1)) - - dzk1=alpha*tanh(l/k0) - dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) - write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2 - zkm1=zk - enddo - - sig(llm+1)=0. - -c - DO 2 l = 1, llm - dsig(l) = sig(l)-sig(l+1) - 2 CONTINUE -c - - ELSE -c----------------------------------------------------------------------- -c cas 2 ancienne discretisation (LMD5...): -c ---------------------------------------- - - WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale' - - if (ok_strato) then - if (llm==39) then - dsigmin=0.3 - else if (llm==50) then - dsigmin=1. - else - WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster' - dsigmin=1. - endif - WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin - endif - - h=7. - snorm = 0. - DO l = 1, llm - x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1) - - IF (ok_strato) THEN - dsig(l) =(dsigmin + 7.0 * SIN(x)**2) - & *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 - ELSE - dsig(l) = 1.0 + 7.0 * SIN(x)**2 - ENDIF - - snorm = snorm + dsig(l) - ENDDO - snorm = 1./snorm - DO l = 1, llm - dsig(l) = dsig(l)*snorm - ENDDO - sig(llm+1) = 0. - DO l = llm, 1, -1 - sig(l) = sig(l+1) + dsig(l) - ENDDO - - ENDIF - - - DO l=1,llm - nivsigs(l) = REAL(l) - ENDDO - - DO l=1,llmp1 - nivsig(l)= REAL(l) - ENDDO - -c -c .... Calculs de ap(l) et de bp(l) .... -c ......................................... -c -c -c ..... pa et preff sont lus sur les fichiers start par lectba ..... -c - - bp(llmp1) = 0. - - DO l = 1, llm -cc -ccc ap(l) = 0. -ccc bp(l) = sig(l) - - bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) - ap(l) = pa * ( sig(l) - bp(l) ) -c - ENDDO - - bp(1)=1. - ap(1)=0. - - ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) - - write(lunout,*)' BP ' - write(lunout,*) bp - write(lunout,*)' AP ' - write(lunout,*) ap - - write(lunout,*) - .'Niveaux de pressions approximatifs aux centres des' - write(lunout,*)'couches calcules pour une pression de surface =', - . preff - write(lunout,*) - . 'et altitudes equivalentes pour une hauteur d echelle de' - write(lunout,*)'8km' - DO l = 1, llm - dpres(l) = bp(l) - bp(l+1) - presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) - write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),' Z ~ ', - . log(preff/presnivs(l))*8. - . ,' DZ ~ ',8.*log((ap(l)+bp(l)*preff)/ - . max(ap(l+1)+bp(l+1)*preff,1.e-10)) - ENDDO - - write(lunout,*)' PRESNIVS ' - write(lunout,*)presnivs - - RETURN - END Index: /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90 =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90 (revision 1472) +++ /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90 (revision 1472) @@ -0,0 +1,142 @@ +! $Id$ + +SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) + + ! Auteur : P. Le Van + + IMPLICIT NONE + + include "dimensions.h" + include "paramet.h" + include "iniprint.h" + include "logic.h" + + ! s = sigma ** kappa : coordonnee verticale + ! dsig(l) : epaisseur de la couche l ds la coord. s + ! sig(l) : sigma a l'interface des couches l et l-1 + ! ds(l) : distance entre les couches l et l-1 en coord.s + + REAL pa, preff + REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) + REAL presnivs(llm) + + REAL sig(llm+1), dsig(llm) + real zk, zkm1, dzk1, dzk2, k0, k1 + + INTEGER l + REAL dsigmin + REAL alpha, beta, deltaz, h + INTEGER iostat + REAL pi, x + + !----------------------------------------------------------------------- + + pi = 2 * ASIN(1.) + + OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) + + IF (iostat == 0) THEN + ! cas 1 on lit les options dans sigma.def: + READ(99, *) h ! hauteur d'echelle 8. + READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 + READ(99, *) beta ! facteur d'acroissement en haut 1.3 + READ(99, *) k0 ! nombre de couches dans la transition surf + READ(99, *) k1 ! nombre de couches dans la transition haute + CLOSE(99) + alpha=deltaz/(llm*h) + write(lunout, *)'h, alpha, k0, k1, beta' + + alpha=deltaz/tanh(1./k0)*2. + zkm1=0. + sig(1)=1. + do l=1, llm + sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) & + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) + zk=-h*log(sig(l+1)) + + dzk1=alpha*tanh(l/k0) + dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) + write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 + zkm1=zk + enddo + + sig(llm+1)=0. + + DO l = 1, llm + dsig(l) = sig(l)-sig(l+1) + end DO + ELSE + if (ok_strato) then + if (llm==39) then + dsigmin=0.3 + else if (llm==50) then + dsigmin=1. + else + WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster' + dsigmin=1. + endif + WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin + endif + + DO l = 1, llm + x = pi * (l - 0.5) / (llm + 1) + + IF (ok_strato) THEN + dsig(l) =(dsigmin + 7. * SIN(x)**2) & + *(1. - tanh(2 * x / pi - 1.))**2 / 4. + ELSE + dsig(l) = 1.0 + 7.0 * SIN(x)**2 + ENDIF + ENDDO + dsig = dsig / sum(dsig) + sig(llm+1) = 0. + DO l = llm, 1, -1 + sig(l) = sig(l+1) + dsig(l) + ENDDO + ENDIF + + DO l=1, llm + nivsigs(l) = REAL(l) + ENDDO + + DO l=1, llmp1 + nivsig(l)= REAL(l) + ENDDO + + ! .... Calculs de ap(l) et de bp(l) .... + ! ..... pa et preff sont lus sur les fichiers start par lectba ..... + + bp(llmp1) = 0. + + DO l = 1, llm + bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) + ap(l) = pa * ( sig(l) - bp(l) ) + ENDDO + + bp(1)=1. + ap(1)=0. + + ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) + + write(lunout, *)' BP ' + write(lunout, *) bp + write(lunout, *)' AP ' + write(lunout, *) ap + + write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' + write(lunout, *)'couches calcules pour une pression de surface =', preff + write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de' + write(lunout, *)'8km' + DO l = 1, llm + dpres(l) = bp(l) - bp(l+1) + presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) + write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & + log(preff/presnivs(l))*8. & + , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ & + max(ap(l+1)+bp(l+1)*preff, 1.e-10)) + ENDDO + + write(lunout, *) 'PRESNIVS ' + write(lunout, *) presnivs + +END SUBROUTINE disvert Index: DZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F (revision 1471) +++ (revision ) @@ -1,285 +1,0 @@ -! -! $Id$ -! -c -c - SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) - - USE filtreg_mod - USE infotrac, ONLY : nqtot - USE control_mod, ONLY: day_step,planet_type -#ifdef CPP_IOIPSL - USE IOIPSL -#else -! if not using IOIPSL, we still need to use (a local version of) getin - USE ioipsl_getincom -#endif - USE Write_Field - - -c%W% %G% -c======================================================================= -c -c Author: Frederic Hourdin original: 15/01/93 -c ------- -c -c Subject: -c ------ -c -c Method: -c -------- -c -c Interface: -c ---------- -c -c Input: -c ------ -c -c Output: -c ------- -c -c======================================================================= - IMPLICIT NONE -c----------------------------------------------------------------------- -c Declararations: -c --------------- - -#include "dimensions.h" -#include "paramet.h" -#include "comvert.h" -#include "comconst.h" -#include "comgeom.h" -#include "academic.h" -#include "ener.h" -#include "temps.h" -#include "iniprint.h" -#include "logic.h" - -c Arguments: -c ---------- - - real time_0 - -c variables dynamiques - REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants - REAL teta(ip1jmp1,llm) ! temperature potentielle - REAL q(ip1jmp1,llm,nqtot) ! champs advectes - REAL ps(ip1jmp1) ! pression au sol - REAL masse(ip1jmp1,llm) ! masse d'air - REAL phis(ip1jmp1) ! geopotentiel au sol - -c Local: -c ------ - - REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches - REAL pks(ip1jmp1) ! exner au sol - REAL pk(ip1jmp1,llm) ! exner au milieu des couches - REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches - REAL phi(ip1jmp1,llm) ! geopotentiel - REAL ddsin,tetastrat,zsig,tetapv,w_pv ! variables auxiliaires - real tetajl(jjp1,llm) - INTEGER i,j,l,lsup,ij - - REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T - REAL k_f,k_c_a,k_c_s ! Constantes de rappel - LOGICAL ok_geost ! Initialisation vent geost. ou nul - LOGICAL ok_pv ! Polar Vortex - REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex - - real zz,ran1 - integer idum - - REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr - -c----------------------------------------------------------------------- -! 1. Initializations for Earth-like case -! -------------------------------------- -c - ! initialize planet radius, rotation rate,... - call conf_planete - - time_0=0. - day_ref=1 - annee_ref=0 - - im = iim - jm = jjm - day_ini = 1 - dtvr = daysec/REAL(day_step) - zdtvr=dtvr - etot0 = 0. - ptot0 = 0. - ztot0 = 0. - stot0 = 0. - ang0 = 0. - - if (llm.eq.1) then - ! specific initializations for the shallow water case - kappa=1 - endif - - CALL iniconst - CALL inigeom - CALL inifilr - - if (llm.eq.1) then - ! initialize fields for the shallow water case, if required - if (.not.read_start) then - phis(:)=0. - q(:,:,:)=0 - CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) - endif - endif - - if (iflag_phys.eq.2) then - ! initializations for the academic case - -! if (planet_type=="earth") then - - ! 1. local parameters - ! by convention, winter is in the southern hemisphere - ! Geostrophic wind or no wind? - ok_geost=.TRUE. - CALL getin('ok_geost',ok_geost) - ! Constants for Newtonian relaxation and friction - k_f=1. !friction - CALL getin('k_j',k_f) - k_f=1./(daysec*k_f) - k_c_s=4. !cooling surface - CALL getin('k_c_s',k_c_s) - k_c_s=1./(daysec*k_c_s) - k_c_a=40. !cooling free atm - CALL getin('k_c_a',k_c_a) - k_c_a=1./(daysec*k_c_a) - ! Constants for Teta equilibrium profile - teta0=315. ! mean Teta (S.H. 315K) - CALL getin('teta0',teta0) - ttp=200. ! Tropopause temperature (S.H. 200K) - CALL getin('ttp',ttp) - eps=0. ! Deviation to N-S symmetry(~0-20K) - CALL getin('eps',eps) - delt_y=60. ! Merid Temp. Gradient (S.H. 60K) - CALL getin('delt_y',delt_y) - delt_z=10. ! Vertical Gradient (S.H. 10K) - CALL getin('delt_z',delt_z) - ! Polar vortex - ok_pv=.false. - CALL getin('ok_pv',ok_pv) - phi_pv=-50. ! Latitude of edge of vortex - CALL getin('phi_pv',phi_pv) - phi_pv=phi_pv*pi/180. - dphi_pv=5. ! Width of the edge - CALL getin('dphi_pv',dphi_pv) - dphi_pv=dphi_pv*pi/180. - gam_pv=4. ! -dT/dz vortex (in K/km) - CALL getin('gam_pv',gam_pv) - - ! 2. Initialize fields towards which to relax - ! Friction - knewt_g=k_c_a - DO l=1,llm - zsig=presnivs(l)/preff - knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) - kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) - ENDDO - DO j=1,jjp1 - clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 - ENDDO - - ! Potential temperature - DO l=1,llm - zsig=presnivs(l)/preff - tetastrat=ttp*zsig**(-kappa) - tetapv=tetastrat - IF ((ok_pv).AND.(zsig.LT.0.1)) THEN - tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) - ENDIF - DO j=1,jjp1 - ! Troposphere - ddsin=sin(rlatu(j)) - tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin - & -delt_z*(1.-ddsin*ddsin)*log(zsig) - ! Profil stratospherique isotherme (+vortex) - w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. - tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv - tetajl(j,l)=MAX(tetajl(j,l),tetastrat) - ENDDO - ENDDO ! of DO l=1,llm - -! CALL writefield('theta_eq',tetajl) - - do l=1,llm - do j=1,jjp1 - do i=1,iip1 - ij=(j-1)*iip1+i - tetarappel(ij,l)=tetajl(j,l) - enddo - enddo - enddo - - -! else -! write(lunout,*)"iniacademic: planet types other than earth", -! & " not implemented (yet)." -! stop -! endif ! of if (planet_type=="earth") - - ! 3. Initialize fields (if necessary) - IF (.NOT. read_start) THEN - ! surface pressure - ps(:)=preff - ! ground geopotential - phis(:)=0. - - CALL pression ( ip1jmp1, ap, bp, ps, p ) - CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) - CALL massdair(p,masse) - - ! bulk initialization of temperature - teta(:,:)=tetarappel(:,:) - - ! geopotential - CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) - - ! winds - if (ok_geost) then - call ugeostr(phi,ucov) - else - ucov(:,:)=0. - endif - vcov(:,:)=0. - - ! bulk initialization of tracers - if (planet_type=="earth") then - ! Earth: first two tracers will be water - do i=1,nqtot - if (i.eq.1) q(:,:,i)=1.e-10 - if (i.eq.2) q(:,:,i)=1.e-15 - if (i.gt.2) q(:,:,i)=0. - enddo - else - q(:,:,:)=0 - endif ! of if (planet_type=="earth") - - ! add random perturbation to temperature - idum = -1 - zz = ran1(idum) - idum = 0 - do l=1,llm - do ij=iip2,ip1jm - teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) - enddo - enddo - - ! maintain periodicity in longitude - do l=1,llm - do ij=1,ip1jmp1,iip1 - teta(ij+iim,l)=teta(ij,l) - enddo - enddo - - ENDIF ! of IF (.NOT. read_start) - endif ! of if (iflag_phys.eq.2) - - END -c----------------------------------------------------------------------- Index: /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90 =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90 (revision 1472) +++ /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90 (revision 1472) @@ -0,0 +1,252 @@ +! +! $Id$ +! +SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) + + USE filtreg_mod + USE infotrac, ONLY : nqtot + USE control_mod, ONLY: day_step,planet_type +#ifdef CPP_IOIPSL + USE IOIPSL +#else + ! if not using IOIPSL, we still need to use (a local version of) getin + USE ioipsl_getincom +#endif + USE Write_Field + + ! Author: Frederic Hourdin original: 15/01/93 + + IMPLICIT NONE + + ! Declararations: + ! --------------- + + include "dimensions.h" + include "paramet.h" + include "comvert.h" + include "comconst.h" + include "comgeom.h" + include "academic.h" + include "ener.h" + include "temps.h" + include "iniprint.h" + include "logic.h" + + ! Arguments: + ! ---------- + + real time_0 + + ! variables dynamiques + REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants + REAL teta(ip1jmp1,llm) ! temperature potentielle + REAL q(ip1jmp1,llm,nqtot) ! champs advectes + REAL ps(ip1jmp1) ! pression au sol + REAL masse(ip1jmp1,llm) ! masse d'air + REAL phis(ip1jmp1) ! geopotentiel au sol + + ! Local: + ! ------ + + REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches + REAL pks(ip1jmp1) ! exner au sol + REAL pk(ip1jmp1,llm) ! exner au milieu des couches + REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches + REAL phi(ip1jmp1,llm) ! geopotentiel + REAL ddsin,tetastrat,zsig,tetapv,w_pv ! variables auxiliaires + real tetajl(jjp1,llm) + INTEGER i,j,l,lsup,ij + + REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T + REAL k_f,k_c_a,k_c_s ! Constantes de rappel + LOGICAL ok_geost ! Initialisation vent geost. ou nul + LOGICAL ok_pv ! Polar Vortex + REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex + + real zz,ran1 + integer idum + + REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr + + !----------------------------------------------------------------------- + ! 1. Initializations for Earth-like case + ! -------------------------------------- + ! + ! initialize planet radius, rotation rate,... + call conf_planete + + time_0=0. + day_ref=1 + annee_ref=0 + + im = iim + jm = jjm + day_ini = 1 + dtvr = daysec/REAL(day_step) + zdtvr=dtvr + etot0 = 0. + ptot0 = 0. + ztot0 = 0. + stot0 = 0. + ang0 = 0. + + if (llm == 1) then + ! specific initializations for the shallow water case + kappa=1 + endif + + CALL iniconst + CALL inigeom + CALL inifilr + + if (llm == 1) then + ! initialize fields for the shallow water case, if required + if (.not.read_start) then + phis(:)=0. + q(:,:,:)=0 + CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) + endif + endif + + academic_case: if (iflag_phys == 2) then + ! initializations + + ! 1. local parameters + ! by convention, winter is in the southern hemisphere + ! Geostrophic wind or no wind? + ok_geost=.TRUE. + CALL getin('ok_geost',ok_geost) + ! Constants for Newtonian relaxation and friction + k_f=1. !friction + CALL getin('k_j',k_f) + k_f=1./(daysec*k_f) + k_c_s=4. !cooling surface + CALL getin('k_c_s',k_c_s) + k_c_s=1./(daysec*k_c_s) + k_c_a=40. !cooling free atm + CALL getin('k_c_a',k_c_a) + k_c_a=1./(daysec*k_c_a) + ! Constants for Teta equilibrium profile + teta0=315. ! mean Teta (S.H. 315K) + CALL getin('teta0',teta0) + ttp=200. ! Tropopause temperature (S.H. 200K) + CALL getin('ttp',ttp) + eps=0. ! Deviation to N-S symmetry(~0-20K) + CALL getin('eps',eps) + delt_y=60. ! Merid Temp. Gradient (S.H. 60K) + CALL getin('delt_y',delt_y) + delt_z=10. ! Vertical Gradient (S.H. 10K) + CALL getin('delt_z',delt_z) + ! Polar vortex + ok_pv=.false. + CALL getin('ok_pv',ok_pv) + phi_pv=-50. ! Latitude of edge of vortex + CALL getin('phi_pv',phi_pv) + phi_pv=phi_pv*pi/180. + dphi_pv=5. ! Width of the edge + CALL getin('dphi_pv',dphi_pv) + dphi_pv=dphi_pv*pi/180. + gam_pv=4. ! -dT/dz vortex (in K/km) + CALL getin('gam_pv',gam_pv) + + ! 2. Initialize fields towards which to relax + ! Friction + knewt_g=k_c_a + DO l=1,llm + zsig=presnivs(l)/preff + knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) + kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) + ENDDO + DO j=1,jjp1 + clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 + ENDDO + + ! Potential temperature + DO l=1,llm + zsig=presnivs(l)/preff + tetastrat=ttp*zsig**(-kappa) + tetapv=tetastrat + IF ((ok_pv).AND.(zsig.LT.0.1)) THEN + tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) + ENDIF + DO j=1,jjp1 + ! Troposphere + ddsin=sin(rlatu(j)) + tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin & + -delt_z*(1.-ddsin*ddsin)*log(zsig) + ! Profil stratospherique isotherme (+vortex) + w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. + tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv + tetajl(j,l)=MAX(tetajl(j,l),tetastrat) + ENDDO + ENDDO + + ! CALL writefield('theta_eq',tetajl) + + do l=1,llm + do j=1,jjp1 + do i=1,iip1 + ij=(j-1)*iip1+i + tetarappel(ij,l)=tetajl(j,l) + enddo + enddo + enddo + + ! 3. Initialize fields (if necessary) + IF (.NOT. read_start) THEN + ! surface pressure + ps(:)=preff + ! ground geopotential + phis(:)=0. + + CALL pression ( ip1jmp1, ap, bp, ps, p ) + CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + CALL massdair(p,masse) + + ! bulk initialization of temperature + teta(:,:)=tetarappel(:,:) + + ! geopotential + CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) + + ! winds + if (ok_geost) then + call ugeostr(phi,ucov) + else + ucov(:,:)=0. + endif + vcov(:,:)=0. + + ! bulk initialization of tracers + if (planet_type=="earth") then + ! Earth: first two tracers will be water + do i=1,nqtot + if (i == 1) q(:,:,i)=1.e-10 + if (i == 2) q(:,:,i)=1.e-15 + if (i.gt.2) q(:,:,i)=0. + enddo + else + q(:,:,:)=0 + endif ! of if (planet_type=="earth") + + ! add random perturbation to temperature + idum = -1 + zz = ran1(idum) + idum = 0 + do l=1,llm + do ij=iip2,ip1jm + teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) + enddo + enddo + + ! maintain periodicity in longitude + do l=1,llm + do ij=1,ip1jmp1,iip1 + teta(ij+iim,l)=teta(ij,l) + enddo + enddo + + ENDIF ! of IF (.NOT. read_start) + endif academic_case + +END SUBROUTINE iniacademic Index: DZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F (revision 1471) +++ (revision ) @@ -1,194 +1,0 @@ -! -! $Id$ -! - SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) - -c Auteur : P. Le Van . -c - IMPLICIT NONE - -#include "dimensions.h" -#include "paramet.h" -#include "iniprint.h" -#include "logic.h" -c -c======================================================================= -c -c -c s = sigma ** kappa : coordonnee verticale -c dsig(l) : epaisseur de la couche l ds la coord. s -c sig(l) : sigma a l'interface des couches l et l-1 -c ds(l) : distance entre les couches l et l-1 en coord.s -c -c======================================================================= -c - REAL pa,preff - REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1) - REAL presnivs(llm) -c -c declarations: -c ------------- -c - REAL sig(llm+1),dsig(llm) - real zzz(1:llm+1) - real dzz(1:llm) - real zk,zkm1,dzk1,dzk2,k0,k1 -c - INTEGER l - REAL snorm,dsigmin - REAL alpha,beta,gama,delta,deltaz,h - INTEGER np,ierr - REAL pi,x - - REAL SSUM -c -c----------------------------------------------------------------------- -c - pi=2.*ASIN(1.) - - OPEN(99,file='sigma.def',status='old',form='formatted', - s iostat=ierr) - -c----------------------------------------------------------------------- -c cas 1 on lit les options dans sigma.def: -c ---------------------------------------- - - IF (ierr.eq.0) THEN - - READ(99,*) h ! hauteur d'echelle 8. - READ(99,*) deltaz ! epaiseur de la premiere couche 0.04 - READ(99,*) beta ! facteur d'acroissement en haut 1.3 - READ(99,*) k0 ! nombre de couches dans la transition surf - READ(99,*) k1 ! nombre de couches dans la transition haute - CLOSE(99) - alpha=deltaz/(llm*h) - write(lunout,*)'h,alpha,k0,k1,beta' - -c read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2 - - alpha=deltaz/tanh(1./k0)*2. - zkm1=0. - sig(1)=1. - do l=1,llm - sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) - + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) - zk=-h*log(sig(l+1)) - - dzk1=alpha*tanh(l/k0) - dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) - write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2 - zkm1=zk - enddo - - sig(llm+1)=0. - -c - DO 2 l = 1, llm - dsig(l) = sig(l)-sig(l+1) - 2 CONTINUE -c - - ELSE -c----------------------------------------------------------------------- -c cas 2 ancienne discretisation (LMD5...): -c ---------------------------------------- - - WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale' - - if (ok_strato) then - if (llm==39) then - dsigmin=0.3 - else if (llm==50) then - dsigmin=1. - else - WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster' - dsigmin=1. - endif - WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin - endif - - h=7. - snorm = 0. - DO l = 1, llm - x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1) - - IF (ok_strato) THEN - dsig(l) =(dsigmin + 7.0 * SIN(x)**2) - & *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 - ELSE - dsig(l) = 1.0 + 7.0 * SIN(x)**2 - ENDIF - - snorm = snorm + dsig(l) - ENDDO - snorm = 1./snorm - DO l = 1, llm - dsig(l) = dsig(l)*snorm - ENDDO - sig(llm+1) = 0. - DO l = llm, 1, -1 - sig(l) = sig(l+1) + dsig(l) - ENDDO - - ENDIF - - - DO l=1,llm - nivsigs(l) = REAL(l) - ENDDO - - DO l=1,llmp1 - nivsig(l)= REAL(l) - ENDDO - -c -c .... Calculs de ap(l) et de bp(l) .... -c ......................................... -c -c -c ..... pa et preff sont lus sur les fichiers start par lectba ..... -c - - bp(llmp1) = 0. - - DO l = 1, llm -cc -ccc ap(l) = 0. -ccc bp(l) = sig(l) - - bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) - ap(l) = pa * ( sig(l) - bp(l) ) -c - ENDDO - - bp(1)=1. - ap(1)=0. - - ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) - - write(lunout,*)' BP ' - write(lunout,*) bp - write(lunout,*)' AP ' - write(lunout,*) ap - - write(lunout,*) - .'Niveaux de pressions approximatifs aux centres des' - write(lunout,*)'couches calcules pour une pression de surface =', - . preff - write(lunout,*) - . 'et altitudes equivalentes pour une hauteur d echelle de' - write(lunout,*)'8km' - DO l = 1, llm - dpres(l) = bp(l) - bp(l+1) - presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) - write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),' Z ~ ', - . log(preff/presnivs(l))*8. - . ,' DZ ~ ',8.*log((ap(l)+bp(l)*preff)/ - . max(ap(l+1)+bp(l+1)*preff,1.e-10)) - ENDDO - - write(lunout,*)' PRESNIVS ' - write(lunout,*)presnivs - - RETURN - END Index: /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90 =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90 (revision 1472) +++ /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90 (revision 1472) @@ -0,0 +1,142 @@ +! $Id$ + +SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig) + + ! Auteur : P. Le Van + + IMPLICIT NONE + + include "dimensions.h" + include "paramet.h" + include "iniprint.h" + include "logic.h" + + ! s = sigma ** kappa : coordonnee verticale + ! dsig(l) : epaisseur de la couche l ds la coord. s + ! sig(l) : sigma a l'interface des couches l et l-1 + ! ds(l) : distance entre les couches l et l-1 en coord.s + + REAL pa, preff + REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1) + REAL presnivs(llm) + + REAL sig(llm+1), dsig(llm) + real zk, zkm1, dzk1, dzk2, k0, k1 + + INTEGER l + REAL dsigmin + REAL alpha, beta, deltaz, h + INTEGER iostat + REAL pi, x + + !----------------------------------------------------------------------- + + pi = 2 * ASIN(1.) + + OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) + + IF (iostat == 0) THEN + ! cas 1 on lit les options dans sigma.def: + READ(99, *) h ! hauteur d'echelle 8. + READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 + READ(99, *) beta ! facteur d'acroissement en haut 1.3 + READ(99, *) k0 ! nombre de couches dans la transition surf + READ(99, *) k1 ! nombre de couches dans la transition haute + CLOSE(99) + alpha=deltaz/(llm*h) + write(lunout, *)'h, alpha, k0, k1, beta' + + alpha=deltaz/tanh(1./k0)*2. + zkm1=0. + sig(1)=1. + do l=1, llm + sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) & + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) + zk=-h*log(sig(l+1)) + + dzk1=alpha*tanh(l/k0) + dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) + write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 + zkm1=zk + enddo + + sig(llm+1)=0. + + DO l = 1, llm + dsig(l) = sig(l)-sig(l+1) + end DO + ELSE + if (ok_strato) then + if (llm==39) then + dsigmin=0.3 + else if (llm==50) then + dsigmin=1. + else + WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster' + dsigmin=1. + endif + WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin + endif + + DO l = 1, llm + x = pi * (l - 0.5) / (llm + 1) + + IF (ok_strato) THEN + dsig(l) =(dsigmin + 7. * SIN(x)**2) & + *(1. - tanh(2 * x / pi - 1.))**2 / 4. + ELSE + dsig(l) = 1.0 + 7.0 * SIN(x)**2 + ENDIF + ENDDO + dsig = dsig / sum(dsig) + sig(llm+1) = 0. + DO l = llm, 1, -1 + sig(l) = sig(l+1) + dsig(l) + ENDDO + ENDIF + + DO l=1, llm + nivsigs(l) = REAL(l) + ENDDO + + DO l=1, llmp1 + nivsig(l)= REAL(l) + ENDDO + + ! .... Calculs de ap(l) et de bp(l) .... + ! ..... pa et preff sont lus sur les fichiers start par lectba ..... + + bp(llmp1) = 0. + + DO l = 1, llm + bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) + ap(l) = pa * ( sig(l) - bp(l) ) + ENDDO + + bp(1)=1. + ap(1)=0. + + ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) + + write(lunout, *)' BP ' + write(lunout, *) bp + write(lunout, *)' AP ' + write(lunout, *) ap + + write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' + write(lunout, *)'couches calcules pour une pression de surface =', preff + write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de' + write(lunout, *)'8km' + DO l = 1, llm + dpres(l) = bp(l) - bp(l+1) + presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) + write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & + log(preff/presnivs(l))*8. & + , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ & + max(ap(l+1)+bp(l+1)*preff, 1.e-10)) + ENDDO + + write(lunout, *) 'PRESNIVS ' + write(lunout, *) presnivs + +END SUBROUTINE disvert Index: DZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F (revision 1471) +++ (revision ) @@ -1,285 +1,0 @@ -! -! $Id$ -! -c -c - SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) - - USE filtreg_mod - USE infotrac, ONLY : nqtot - USE control_mod, ONLY: day_step,planet_type -#ifdef CPP_IOIPSL - USE IOIPSL -#else -! if not using IOIPSL, we still need to use (a local version of) getin - USE ioipsl_getincom -#endif - USE Write_Field - - -c%W% %G% -c======================================================================= -c -c Author: Frederic Hourdin original: 15/01/93 -c ------- -c -c Subject: -c ------ -c -c Method: -c -------- -c -c Interface: -c ---------- -c -c Input: -c ------ -c -c Output: -c ------- -c -c======================================================================= - IMPLICIT NONE -c----------------------------------------------------------------------- -c Declararations: -c --------------- - -#include "dimensions.h" -#include "paramet.h" -#include "comvert.h" -#include "comconst.h" -#include "comgeom.h" -#include "academic.h" -#include "ener.h" -#include "temps.h" -#include "iniprint.h" -#include "logic.h" - -c Arguments: -c ---------- - - real time_0 - -c variables dynamiques - REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants - REAL teta(ip1jmp1,llm) ! temperature potentielle - REAL q(ip1jmp1,llm,nqtot) ! champs advectes - REAL ps(ip1jmp1) ! pression au sol - REAL masse(ip1jmp1,llm) ! masse d'air - REAL phis(ip1jmp1) ! geopotentiel au sol - -c Local: -c ------ - - REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches - REAL pks(ip1jmp1) ! exner au sol - REAL pk(ip1jmp1,llm) ! exner au milieu des couches - REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches - REAL phi(ip1jmp1,llm) ! geopotentiel - REAL ddsin,tetastrat,zsig,tetapv,w_pv ! variables auxiliaires - real tetajl(jjp1,llm) - INTEGER i,j,l,lsup,ij - - REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T - REAL k_f,k_c_a,k_c_s ! Constantes de rappel - LOGICAL ok_geost ! Initialisation vent geost. ou nul - LOGICAL ok_pv ! Polar Vortex - REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex - - real zz,ran1 - integer idum - - REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr - -c----------------------------------------------------------------------- -! 1. Initializations for Earth-like case -! -------------------------------------- -c - ! initialize planet radius, rotation rate,... - call conf_planete - - time_0=0. - day_ref=1 - annee_ref=0 - - im = iim - jm = jjm - day_ini = 1 - dtvr = daysec/REAL(day_step) - zdtvr=dtvr - etot0 = 0. - ptot0 = 0. - ztot0 = 0. - stot0 = 0. - ang0 = 0. - - if (llm.eq.1) then - ! specific initializations for the shallow water case - kappa=1 - endif - - CALL iniconst - CALL inigeom - CALL inifilr - - if (llm.eq.1) then - ! initialize fields for the shallow water case, if required - if (.not.read_start) then - phis(:)=0. - q(:,:,:)=0 - CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) - endif - endif - - if (iflag_phys.eq.2) then - ! initializations for the academic case - -! if (planet_type=="earth") then - - ! 1. local parameters - ! by convention, winter is in the southern hemisphere - ! Geostrophic wind or no wind? - ok_geost=.TRUE. - CALL getin('ok_geost',ok_geost) - ! Constants for Newtonian relaxation and friction - k_f=1. !friction - CALL getin('k_j',k_f) - k_f=1./(daysec*k_f) - k_c_s=4. !cooling surface - CALL getin('k_c_s',k_c_s) - k_c_s=1./(daysec*k_c_s) - k_c_a=40. !cooling free atm - CALL getin('k_c_a',k_c_a) - k_c_a=1./(daysec*k_c_a) - ! Constants for Teta equilibrium profile - teta0=315. ! mean Teta (S.H. 315K) - CALL getin('teta0',teta0) - ttp=200. ! Tropopause temperature (S.H. 200K) - CALL getin('ttp',ttp) - eps=0. ! Deviation to N-S symmetry(~0-20K) - CALL getin('eps',eps) - delt_y=60. ! Merid Temp. Gradient (S.H. 60K) - CALL getin('delt_y',delt_y) - delt_z=10. ! Vertical Gradient (S.H. 10K) - CALL getin('delt_z',delt_z) - ! Polar vortex - ok_pv=.false. - CALL getin('ok_pv',ok_pv) - phi_pv=-50. ! Latitude of edge of vortex - CALL getin('phi_pv',phi_pv) - phi_pv=phi_pv*pi/180. - dphi_pv=5. ! Width of the edge - CALL getin('dphi_pv',dphi_pv) - dphi_pv=dphi_pv*pi/180. - gam_pv=4. ! -dT/dz vortex (in K/km) - CALL getin('gam_pv',gam_pv) - - ! 2. Initialize fields towards which to relax - ! Friction - knewt_g=k_c_a - DO l=1,llm - zsig=presnivs(l)/preff - knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) - kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) - ENDDO - DO j=1,jjp1 - clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 - ENDDO - - ! Potential temperature - DO l=1,llm - zsig=presnivs(l)/preff - tetastrat=ttp*zsig**(-kappa) - tetapv=tetastrat - IF ((ok_pv).AND.(zsig.LT.0.1)) THEN - tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) - ENDIF - DO j=1,jjp1 - ! Troposphere - ddsin=sin(rlatu(j)) - tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin - & -delt_z*(1.-ddsin*ddsin)*log(zsig) - ! Profil stratospherique isotherme (+vortex) - w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. - tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv - tetajl(j,l)=MAX(tetajl(j,l),tetastrat) - ENDDO - ENDDO ! of DO l=1,llm - -! CALL writefield('theta_eq',tetajl) - - do l=1,llm - do j=1,jjp1 - do i=1,iip1 - ij=(j-1)*iip1+i - tetarappel(ij,l)=tetajl(j,l) - enddo - enddo - enddo - - -! else -! write(lunout,*)"iniacademic: planet types other than earth", -! & " not implemented (yet)." -! stop -! endif ! of if (planet_type=="earth") - - ! 3. Initialize fields (if necessary) - IF (.NOT. read_start) THEN - ! surface pressure - ps(:)=preff - ! ground geopotential - phis(:)=0. - - CALL pression ( ip1jmp1, ap, bp, ps, p ) - CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) - CALL massdair(p,masse) - - ! bulk initialization of temperature - teta(:,:)=tetarappel(:,:) - - ! geopotential - CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) - - ! winds - if (ok_geost) then - call ugeostr(phi,ucov) - else - ucov(:,:)=0. - endif - vcov(:,:)=0. - - ! bulk initialization of tracers - if (planet_type=="earth") then - ! Earth: first two tracers will be water - do i=1,nqtot - if (i.eq.1) q(:,:,i)=1.e-10 - if (i.eq.2) q(:,:,i)=1.e-15 - if (i.gt.2) q(:,:,i)=0. - enddo - else - q(:,:,:)=0 - endif ! of if (planet_type=="earth") - - ! add random perturbation to temperature - idum = -1 - zz = ran1(idum) - idum = 0 - do l=1,llm - do ij=iip2,ip1jm - teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) - enddo - enddo - - ! maintain periodicity in longitude - do l=1,llm - do ij=1,ip1jmp1,iip1 - teta(ij+iim,l)=teta(ij,l) - enddo - enddo - - ENDIF ! of IF (.NOT. read_start) - endif ! of if (iflag_phys.eq.2) - - END -c----------------------------------------------------------------------- Index: /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F90 =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F90 (revision 1472) +++ /LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F90 (revision 1472) @@ -0,0 +1,252 @@ +! +! $Id$ +! +SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) + + USE filtreg_mod + USE infotrac, ONLY : nqtot + USE control_mod, ONLY: day_step,planet_type +#ifdef CPP_IOIPSL + USE IOIPSL +#else + ! if not using IOIPSL, we still need to use (a local version of) getin + USE ioipsl_getincom +#endif + USE Write_Field + + ! Author: Frederic Hourdin original: 15/01/93 + + IMPLICIT NONE + + ! Declararations: + ! --------------- + + include "dimensions.h" + include "paramet.h" + include "comvert.h" + include "comconst.h" + include "comgeom.h" + include "academic.h" + include "ener.h" + include "temps.h" + include "iniprint.h" + include "logic.h" + + ! Arguments: + ! ---------- + + real time_0 + + ! variables dynamiques + REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants + REAL teta(ip1jmp1,llm) ! temperature potentielle + REAL q(ip1jmp1,llm,nqtot) ! champs advectes + REAL ps(ip1jmp1) ! pression au sol + REAL masse(ip1jmp1,llm) ! masse d'air + REAL phis(ip1jmp1) ! geopotentiel au sol + + ! Local: + ! ------ + + REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches + REAL pks(ip1jmp1) ! exner au sol + REAL pk(ip1jmp1,llm) ! exner au milieu des couches + REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches + REAL phi(ip1jmp1,llm) ! geopotentiel + REAL ddsin,tetastrat,zsig,tetapv,w_pv ! variables auxiliaires + real tetajl(jjp1,llm) + INTEGER i,j,l,lsup,ij + + REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T + REAL k_f,k_c_a,k_c_s ! Constantes de rappel + LOGICAL ok_geost ! Initialisation vent geost. ou nul + LOGICAL ok_pv ! Polar Vortex + REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex + + real zz,ran1 + integer idum + + REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr + + !----------------------------------------------------------------------- + ! 1. Initializations for Earth-like case + ! -------------------------------------- + ! + ! initialize planet radius, rotation rate,... + call conf_planete + + time_0=0. + day_ref=1 + annee_ref=0 + + im = iim + jm = jjm + day_ini = 1 + dtvr = daysec/REAL(day_step) + zdtvr=dtvr + etot0 = 0. + ptot0 = 0. + ztot0 = 0. + stot0 = 0. + ang0 = 0. + + if (llm == 1) then + ! specific initializations for the shallow water case + kappa=1 + endif + + CALL iniconst + CALL inigeom + CALL inifilr + + if (llm == 1) then + ! initialize fields for the shallow water case, if required + if (.not.read_start) then + phis(:)=0. + q(:,:,:)=0 + CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) + endif + endif + + academic_case: if (iflag_phys == 2) then + ! initializations + + ! 1. local parameters + ! by convention, winter is in the southern hemisphere + ! Geostrophic wind or no wind? + ok_geost=.TRUE. + CALL getin('ok_geost',ok_geost) + ! Constants for Newtonian relaxation and friction + k_f=1. !friction + CALL getin('k_j',k_f) + k_f=1./(daysec*k_f) + k_c_s=4. !cooling surface + CALL getin('k_c_s',k_c_s) + k_c_s=1./(daysec*k_c_s) + k_c_a=40. !cooling free atm + CALL getin('k_c_a',k_c_a) + k_c_a=1./(daysec*k_c_a) + ! Constants for Teta equilibrium profile + teta0=315. ! mean Teta (S.H. 315K) + CALL getin('teta0',teta0) + ttp=200. ! Tropopause temperature (S.H. 200K) + CALL getin('ttp',ttp) + eps=0. ! Deviation to N-S symmetry(~0-20K) + CALL getin('eps',eps) + delt_y=60. ! Merid Temp. Gradient (S.H. 60K) + CALL getin('delt_y',delt_y) + delt_z=10. ! Vertical Gradient (S.H. 10K) + CALL getin('delt_z',delt_z) + ! Polar vortex + ok_pv=.false. + CALL getin('ok_pv',ok_pv) + phi_pv=-50. ! Latitude of edge of vortex + CALL getin('phi_pv',phi_pv) + phi_pv=phi_pv*pi/180. + dphi_pv=5. ! Width of the edge + CALL getin('dphi_pv',dphi_pv) + dphi_pv=dphi_pv*pi/180. + gam_pv=4. ! -dT/dz vortex (in K/km) + CALL getin('gam_pv',gam_pv) + + ! 2. Initialize fields towards which to relax + ! Friction + knewt_g=k_c_a + DO l=1,llm + zsig=presnivs(l)/preff + knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) + kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) + ENDDO + DO j=1,jjp1 + clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 + ENDDO + + ! Potential temperature + DO l=1,llm + zsig=presnivs(l)/preff + tetastrat=ttp*zsig**(-kappa) + tetapv=tetastrat + IF ((ok_pv).AND.(zsig.LT.0.1)) THEN + tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) + ENDIF + DO j=1,jjp1 + ! Troposphere + ddsin=sin(rlatu(j)) + tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin & + -delt_z*(1.-ddsin*ddsin)*log(zsig) + ! Profil stratospherique isotherme (+vortex) + w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. + tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv + tetajl(j,l)=MAX(tetajl(j,l),tetastrat) + ENDDO + ENDDO + + ! CALL writefield('theta_eq',tetajl) + + do l=1,llm + do j=1,jjp1 + do i=1,iip1 + ij=(j-1)*iip1+i + tetarappel(ij,l)=tetajl(j,l) + enddo + enddo + enddo + + ! 3. Initialize fields (if necessary) + IF (.NOT. read_start) THEN + ! surface pressure + ps(:)=preff + ! ground geopotential + phis(:)=0. + + CALL pression ( ip1jmp1, ap, bp, ps, p ) + CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) + CALL massdair(p,masse) + + ! bulk initialization of temperature + teta(:,:)=tetarappel(:,:) + + ! geopotential + CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) + + ! winds + if (ok_geost) then + call ugeostr(phi,ucov) + else + ucov(:,:)=0. + endif + vcov(:,:)=0. + + ! bulk initialization of tracers + if (planet_type=="earth") then + ! Earth: first two tracers will be water + do i=1,nqtot + if (i == 1) q(:,:,i)=1.e-10 + if (i == 2) q(:,:,i)=1.e-15 + if (i.gt.2) q(:,:,i)=0. + enddo + else + q(:,:,:)=0 + endif ! of if (planet_type=="earth") + + ! add random perturbation to temperature + idum = -1 + zz = ran1(idum) + idum = 0 + do l=1,llm + do ij=iip2,ip1jm + teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) + enddo + enddo + + ! maintain periodicity in longitude + do l=1,llm + do ij=1,ip1jmp1,iip1 + teta(ij+iim,l)=teta(ij,l) + enddo + enddo + + ENDIF ! of IF (.NOT. read_start) + endif academic_case + +END SUBROUTINE iniacademic Index: DZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F (revision 1471) +++ (revision ) @@ -1,612 +1,0 @@ -! -! $Id$ -! -c - SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, - s d_t, d_q, d_ql, rneb, radliq, rain, snow, - s pfrac_impa, pfrac_nucl, pfrac_1nucl, - s frac_impa, frac_nucl, - s prfl, psfl, rhcl, zqta, fraca, - s ztv, zpspsk, ztla, zthl, iflag_cldcon) - -c - USE dimphy - IMPLICIT none -c====================================================================== -c Auteur(s): Z.X. Li (LMD/CNRS) -c Date: le 20 mars 1995 -c Objet: condensation et precipitation stratiforme. -c schema de nuage -c====================================================================== -c====================================================================== -cym#include "dimensions.h" -cym#include "dimphy.h" -#include "YOMCST.h" -#include "tracstoke.h" -#include "fisrtilp.h" -c -c Arguments: -c - REAL dtime ! intervalle du temps (s) - REAL paprs(klon,klev+1) ! pression a inter-couche - REAL pplay(klon,klev) ! pression au milieu de couche - REAL t(klon,klev) ! temperature (K) - REAL q(klon,klev) ! humidite specifique (kg/kg) - REAL d_t(klon,klev) ! incrementation de la temperature (K) - REAL d_q(klon,klev) ! incrementation de la vapeur d'eau - REAL d_ql(klon,klev) ! incrementation de l'eau liquide - REAL rneb(klon,klev) ! fraction nuageuse - REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements - REAL rhcl(klon,klev) ! humidite relative en ciel clair - REAL rain(klon) ! pluies (mm/s) - REAL snow(klon) ! neige (mm/s) - REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) - REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) - REAL ztv(klon,klev) - REAL zqta(klon,klev),fraca(klon,klev) - REAL sigma1(klon,klev),sigma2(klon,klev) - REAL qltot(klon,klev),ctot(klon,klev) - REAL zpspsk(klon,klev),ztla(klon,klev) - REAL zthl(klon,klev) - - logical lognormale(klon) - -cAA -c Coeffients de fraction lessivee : pour OFF-LINE -c - REAL pfrac_nucl(klon,klev) - REAL pfrac_1nucl(klon,klev) - REAL pfrac_impa(klon,klev) -c -c Fraction d'aerosols lessivee par impaction et par nucleation -c POur ON-LINE -c - REAL frac_impa(klon,klev) - REAL frac_nucl(klon,klev) - real zct ,zcl -cAA -c -c Options du programme: -c - REAL seuil_neb ! un nuage existe vraiment au-dela - PARAMETER (seuil_neb=0.001) - - INTEGER ninter ! sous-intervals pour la precipitation - INTEGER ncoreczq - INTEGER iflag_cldcon - PARAMETER (ninter=5) - LOGICAL evap_prec ! evaporation de la pluie - PARAMETER (evap_prec=.TRUE.) - REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur - logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur - - real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) - real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) - real erf - REAL qcloud(klon) -c - LOGICAL cpartiel ! condensation partielle - PARAMETER (cpartiel=.TRUE.) - REAL t_coup - PARAMETER (t_coup=234.0) -c -c Variables locales: -c - INTEGER i, k, n, kk - REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 - REAL zrfl(klon), zrfln(klon), zqev, zqevt - REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq - REAL ztglace, zt(klon) - INTEGER nexpo ! exponentiel pour glace/eau - REAL zdz(klon),zrho(klon),ztot , zrhol(klon) - REAL zchau ,zfroi ,zfice(klon),zneb(klon) -c - LOGICAL appel1er - SAVE appel1er -c$OMP THREADPRIVATE(appel1er) -c -c--------------------------------------------------------------- -c -cAA Variables traceurs: -cAA Provisoire !!! Parametres alpha du lessivage -cAA A priori on a 4 scavenging # possibles -c - REAL a_tr_sca(4) - save a_tr_sca -c$OMP THREADPRIVATE(a_tr_sca) -c -c Variables intermediaires -c - REAL zalpha_tr - REAL zfrac_lessi - REAL zprec_cond(klon) -cAA - REAL zmair, zcpair, zcpeau -C Pour la conversion eau-neige - REAL zlh_solid(klon), zm_solid -cIM -cym INTEGER klevm1 -c--------------------------------------------------------------- -c -c Fonctions en ligne: -c - REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace - REAL zzz -#include "YOETHF.h" -#include "FCTTRE.h" - fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con - fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc -c - DATA appel1er /.TRUE./ -cym - zdelq=0.0 - - print*,'NUAGES4 A. JAM' - IF (appel1er) THEN -c - PRINT*, 'fisrtilp, ninter:', ninter - PRINT*, 'fisrtilp, evap_prec:', evap_prec - PRINT*, 'fisrtilp, cpartiel:', cpartiel - IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN - PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime - PRINT*, 'Je prefere un sous-intervalle de 6 minutes' -c CALL abort - ENDIF - appel1er = .FALSE. -c -cAA initialiation provisoire - a_tr_sca(1) = -0.5 - a_tr_sca(2) = -0.5 - a_tr_sca(3) = -0.5 - a_tr_sca(4) = -0.5 -c -cAA Initialisation a 1 des coefs des fractions lessivees -c -!cdir collapse - DO k = 1, klev - DO i = 1, klon - pfrac_nucl(i,k)=1. - pfrac_1nucl(i,k)=1. - pfrac_impa(i,k)=1. - ENDDO - ENDDO - - ENDIF ! test sur appel1er -c -cMAf Initialisation a 0 de zoliq -c DO i = 1, klon -c zoliq(i)=0. -c ENDDO -c Determiner les nuages froids par leur temperature -c nexpo regle la raideur de la transition eau liquide / eau glace. -c - ztglace = RTT - 15.0 - nexpo = 6 -ccc nexpo = 1 -c -c Initialiser les sorties: -c -!cdir collapse - DO k = 1, klev+1 - DO i = 1, klon - prfl(i,k) = 0.0 - psfl(i,k) = 0.0 - ENDDO - ENDDO - -!cdir collapse - DO k = 1, klev - DO i = 1, klon - d_t(i,k) = 0.0 - d_q(i,k) = 0.0 - d_ql(i,k) = 0.0 - rneb(i,k) = 0.0 - radliq(i,k) = 0.0 - frac_nucl(i,k) = 1. - frac_impa(i,k) = 1. - ENDDO - ENDDO - DO i = 1, klon - rain(i) = 0.0 - snow(i) = 0.0 - zoliq(i)=0. -c ENDDO -c -c Initialiser le flux de precipitation a zero -c -c DO i = 1, klon - zrfl(i) = 0.0 - zneb(i) = seuil_neb - ENDDO -c -c -cAA Pour plus de securite - - zalpha_tr = 0. - zfrac_lessi = 0. - -cAA---------------------------------------------------------- -c - ncoreczq=0 -c Boucle verticale (du haut vers le bas) -c -cIM : klevm1 -cym klevm1=klev-1 - DO 9999 k = klev, 1, -1 -c -cAA---------------------------------------------------------- -c - DO i = 1, klon - zt(i)=t(i,k) - zq(i)=q(i,k) - ENDDO -c -c Calculer la varition de temp. de l'air du a la chaleur sensible -C transporter par la pluie. -C Il resterait a rajouter cet effet de la chaleur sensible sur les -C flux de surface, du a la diff. de temp. entre le 1er niveau et la -C surface. -C - IF(k.LE.klevm1) THEN - DO i = 1, klon -cIM - zmair=(paprs(i,k)-paprs(i,k+1))/RG - zcpair=RCPD*(1.0+RVTMP2*zq(i)) - zcpeau=RCPD*RVTMP2 - zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau - $ + zmair*zcpair*zt(i) ) - $ / (zmair*zcpair + zrfl(i)*dtime*zcpeau) -C C WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1)) - ENDDO - ENDIF -c -c -c Calculer l'evaporation de la precipitation -c - - - IF (evap_prec) THEN - DO i = 1, klon - IF (zrfl(i) .GT.0.) THEN - IF (thermcep) THEN - zdelta=MAX(0.,SIGN(1.,RTT-zt(i))) - zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k) - zqs(i)=MIN(0.5,zqs(i)) - zcor=1./(1.-RETV*zqs(i)) - zqs(i)=zqs(i)*zcor - ELSE - IF (zt(i) .LT. t_coup) THEN - zqs(i) = qsats(zt(i)) / pplay(i,k) - ELSE - zqs(i) = qsatl(zt(i)) / pplay(i,k) - ENDIF - ENDIF - zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) ) - zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) - . * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG - zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) - . * RG*dtime/(paprs(i,k)-paprs(i,k+1)) - zqev = MIN (zqev, zqevt) - zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) - . /RG/dtime - -c pour la glace, on r�vapore toute la pr�ip dans la couche du dessous -c la glace venant de la couche du dessus est simplement dans la couche -c du dessous. - - IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0. - - zq(i) = zq(i) - (zrfln(i)-zrfl(i)) - . * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime - zt(i) = zt(i) + (zrfln(i)-zrfl(i)) - . * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime - . * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) - zrfl(i) = zrfln(i) - ENDIF - ENDDO - ENDIF -c -c Calculer Qs et L/Cp*dQs/dT: -c - IF (thermcep) THEN - DO i = 1, klon - zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) - zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta - zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i)) - zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k) - zqs(i) = MIN(0.5,zqs(i)) - zcor = 1./(1.-RETV*zqs(i)) - zqs(i) = zqs(i)*zcor - zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor) - ENDDO - ELSE - DO i = 1, klon - IF (zt(i).LT.t_coup) THEN - zqs(i) = qsats(zt(i))/pplay(i,k) - zdqs(i) = dqsats(zt(i),zqs(i)) - ELSE - zqs(i) = qsatl(zt(i))/pplay(i,k) - zdqs(i) = dqsatl(zt(i),zqs(i)) - ENDIF - ENDDO - ENDIF -c -c Determiner la condensation partielle et calculer la quantite -c de l'eau condensee: -c - - IF (cpartiel) THEN - -c print*,'Dans partiel k=',k -c -c Calcul de l'eau condensee et de la fraction nuageuse et de l'eau -c nuageuse a partir des PDF de Sandrine Bony. -c rneb : fraction nuageuse -c zqn : eau totale dans le nuage -c zcond : eau condensee moyenne dans la maille. -c on prend en compte le r�hauffement qui diminue la partie condensee -c -c Version avec les raqts - - if (iflag_pdf.eq.0) then - - do i=1,klon - zdelq = min(ratqs(i,k),0.99) * zq(i) - rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq) - zqn(i) = (zq(i)+zdelq+zqs(i))/2.0 - enddo - - else -c -c Version avec les nouvelles PDFs. - do i=1,klon - if(zq(i).lt.1.e-15) then - ncoreczq=ncoreczq+1 - zq(i)=1.e-15 - endif - enddo - - if (iflag_cldcon>=5) then - - call cloudth(klon,klev,k,ztv, - . zq,zqta,fraca, - . qcloud,ctot,zpspsk,paprs,ztla,zthl, - . ratqs,zqs,t) - - do i=1,klon - rneb(i,k)=ctot(i,k) - zqn(i)=qcloud(i) - enddo - - endif - -! Pour iflag_cldcon<=4, on prend toujours la lognormale -! Dans le cas iflag_cldcon=5, on prend systématiquement la bi-gaussienne -! Dans le cas iflagÃ_cldcon=6, on prend la lognormale en absence des thermiques - - lognormale(:)= - . iflag_cldcon<=4.or.(iflag_cldcon==6.and.fraca(:,k)<1.e-10) - do i=1,klon - if (lognormale(i)) then - zpdf_sig(i)=ratqs(i,k)*zq(i) - zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2)) - zpdf_delta(i)=log(zq(i)/zqs(i)) - zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.)) - zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.)) - zpdf_e1(i)=zpdf_a(i)-zpdf_b(i) - zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i)) - zpdf_e1(i)=1.-erf(zpdf_e1(i)) - zpdf_e2(i)=zpdf_a(i)+zpdf_b(i) - zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i)) - zpdf_e2(i)=1.-erf(zpdf_e2(i)) - endif - enddo - - do i=1,klon - if (lognormale(i)) then - if (zpdf_e1(i).lt.1.e-10) then - rneb(i,k)=0. - zqn(i)=zqs(i) - else - rneb(i,k)=0.5*zpdf_e1(i) - zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i) - endif - endif - - enddo - - - endif ! iflag_pdf - - DO i=1,klon - IF (rneb(i,k) .LE. 0.0) THEN - zqn(i) = 0.0 - rneb(i,k) = 0.0 - zcond(i) = 0.0 - rhcl(i,k)=zq(i)/zqs(i) - ELSE IF (rneb(i,k) .GE. 1.0) THEN - zqn(i) = zq(i) - rneb(i,k) = 1.0 - zcond(i) = MAX(0.0,zqn(i)-zqs(i)) - rhcl(i,k)=1.0 - ELSE - zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k) - rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i) - ENDIF - ENDDO -! do i=1,klon -! IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0 -! IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i) -! rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k))) -!c zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i)) -!c On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par -!c la convection. -!c ATTENTION !!! Il va falloir verifier tout ca. -! zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k) -!c print*,'ZDQS ',zdqs(i) -!c--Olivier -! rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i) -! IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i) -! IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0 -!c--fin -! ENDDO - ELSE - DO i = 1, klon - IF (zq(i).GT.zqs(i)) THEN - rneb(i,k) = 1.0 - ELSE - rneb(i,k) = 0.0 - ENDIF - zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i)) - ENDDO - ENDIF -c - DO i = 1, klon - zq(i) = zq(i) - zcond(i) -c zt(i) = zt(i) + zcond(i) * RLVTT/RCPD - zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) - ENDDO -c -c Partager l'eau condensee en precipitation et eau liquide nuageuse -c - DO i = 1, klon - IF (rneb(i,k).GT.0.0) THEN - zoliq(i) = zcond(i) - zrho(i) = pplay(i,k) / zt(i) / RD - zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG) - zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace) - zfice(i) = MIN(MAX(zfice(i),0.0),1.0) - zfice(i) = zfice(i)**nexpo - zneb(i) = MAX(rneb(i,k), seuil_neb) - radliq(i,k) = zoliq(i)/REAL(ninter+1) - ENDIF - ENDDO -c - DO n = 1, ninter - DO i = 1, klon - IF (rneb(i,k).GT.0.0) THEN - zrhol(i) = zrho(i) * zoliq(i) / zneb(i) - - IF (zneb(i).EQ.seuil_neb) THEN - ztot = 0.0 - ELSE -c quantite d'eau a eliminer: zchau -c meme chose pour la glace: zfroi - if (ptconv(i,k)) then - zcl =cld_lc_con - zct =1./cld_tau_con - zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) - . *fallvc(zrhol(i)) * zfice(i) - else - zcl =cld_lc_lsc - zct =1./cld_tau_lsc - zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) - . *fallvs(zrhol(i)) * zfice(i) - endif - zchau = zct *dtime/REAL(ninter) * zoliq(i) - . *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i)) - ztot = zchau + zfroi - ztot = MAX(ztot ,0.0) - ENDIF - ztot = MIN(ztot,zoliq(i)) - zoliq(i) = MAX(zoliq(i)-ztot , 0.0) - radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1) - ENDIF - ENDDO - ENDDO -c - DO i = 1, klon - IF (rneb(i,k).GT.0.0) THEN - d_ql(i,k) = zoliq(i) - zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) - . * (paprs(i,k)-paprs(i,k+1))/(RG*dtime) - ENDIF - IF (zt(i).LT.RTT) THEN - psfl(i,k)=zrfl(i) - ELSE - prfl(i,k)=zrfl(i) - ENDIF - ENDDO -c -c Calculer les tendances de q et de t: -c - DO i = 1, klon - d_q(i,k) = zq(i) - q(i,k) - d_t(i,k) = zt(i) - t(i,k) - ENDDO -c -cAA--------------- Calcul du lessivage stratiforme ------------- - - DO i = 1,klon -c - zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) - . * (paprs(i,k)-paprs(i,k+1))/RG - IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN -cAA lessivage nucleation LMD5 dans la couche elle-meme - if (t(i,k) .GE. ztglace) THEN - zalpha_tr = a_tr_sca(3) - else - zalpha_tr = a_tr_sca(4) - endif - zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) - pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi) - frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi -c -c nucleation avec un facteur -1 au lieu de -0.5 - zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i)) - pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi) - ENDIF -c - ENDDO ! boucle sur i -c -cAA Lessivage par impaction dans les couches en-dessous - DO kk = k-1, 1, -1 - DO i = 1, klon - IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN - if (t(i,kk) .GE. ztglace) THEN - zalpha_tr = a_tr_sca(1) - else - zalpha_tr = a_tr_sca(2) - endif - zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) - pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi) - frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi - ENDIF - ENDDO - ENDDO -c -cAA---------------------------------------------------------- -c FIN DE BOUCLE SUR K - 9999 CONTINUE -c -cAA----------------------------------------------------------- -c -c Pluie ou neige au sol selon la temperature de la 1ere couche -c - DO i = 1, klon - IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN - snow(i) = zrfl(i) - zlh_solid(i) = RLSTT-RLVTT - ELSE - rain(i) = zrfl(i) - zlh_solid(i) = 0. - ENDIF - ENDDO -C -C For energy conservation : when snow is present, the solification -c latent heat is considered. - DO k = 1, klev - DO i = 1, klon - zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k))) - zmair=(paprs(i,k)-paprs(i,k+1))/RG - zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime - d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair) - END DO - END DO -c - - if (ncoreczq>0) then - print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.' - endif - RETURN - END Index: /LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F90 =================================================================== --- /LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F90 (revision 1472) +++ /LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F90 (revision 1472) @@ -0,0 +1,619 @@ +! +! $Id$ +! +! +SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, & + d_t, d_q, d_ql, rneb, radliq, rain, snow, & + pfrac_impa, pfrac_nucl, pfrac_1nucl, & + frac_impa, frac_nucl, & + prfl, psfl, rhcl, zqta, fraca, & + ztv, zpspsk, ztla, zthl, iflag_cldcon) + + ! + USE dimphy + IMPLICIT none + !====================================================================== + ! Auteur(s): Z.X. Li (LMD/CNRS) + ! Date: le 20 mars 1995 + ! Objet: condensation et precipitation stratiforme. + ! schema de nuage + !====================================================================== + !====================================================================== + !ym include "dimensions.h" + !ym include "dimphy.h" + include "YOMCST.h" + include "tracstoke.h" + include "fisrtilp.h" + ! + ! Arguments: + ! + REAL dtime ! intervalle du temps (s) + REAL paprs(klon,klev+1) ! pression a inter-couche + REAL pplay(klon,klev) ! pression au milieu de couche + REAL t(klon,klev) ! temperature (K) + REAL q(klon,klev) ! humidite specifique (kg/kg) + REAL d_t(klon,klev) ! incrementation de la temperature (K) + REAL d_q(klon,klev) ! incrementation de la vapeur d'eau + REAL d_ql(klon,klev) ! incrementation de l'eau liquide + REAL rneb(klon,klev) ! fraction nuageuse + REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements + REAL rhcl(klon,klev) ! humidite relative en ciel clair + REAL rain(klon) ! pluies (mm/s) + REAL snow(klon) ! neige (mm/s) + REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) + REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) + REAL ztv(klon,klev) + REAL zqta(klon,klev),fraca(klon,klev) + REAL sigma1(klon,klev),sigma2(klon,klev) + REAL qltot(klon,klev),ctot(klon,klev) + REAL zpspsk(klon,klev),ztla(klon,klev) + REAL zthl(klon,klev) + + logical lognormale(klon) + + !AA + ! Coeffients de fraction lessivee : pour OFF-LINE + ! + REAL pfrac_nucl(klon,klev) + REAL pfrac_1nucl(klon,klev) + REAL pfrac_impa(klon,klev) + ! + ! Fraction d'aerosols lessivee par impaction et par nucleation + ! POur ON-LINE + ! + REAL frac_impa(klon,klev) + REAL frac_nucl(klon,klev) + real zct ,zcl + !AA + ! + ! Options du programme: + ! + REAL seuil_neb ! un nuage existe vraiment au-dela + PARAMETER (seuil_neb=0.001) + + INTEGER ninter ! sous-intervals pour la precipitation + INTEGER ncoreczq + INTEGER iflag_cldcon + PARAMETER (ninter=5) + LOGICAL evap_prec ! evaporation de la pluie + PARAMETER (evap_prec=.TRUE.) + REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur + logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur + + real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) + real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) + real erf + REAL qcloud(klon) + ! + LOGICAL cpartiel ! condensation partielle + PARAMETER (cpartiel=.TRUE.) + REAL t_coup + PARAMETER (t_coup=234.0) + ! + ! Variables locales: + ! + INTEGER i, k, n, kk + REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 + REAL zrfl(klon), zrfln(klon), zqev, zqevt + REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq + REAL ztglace, zt(klon) + INTEGER nexpo ! exponentiel pour glace/eau + REAL zdz(klon),zrho(klon),ztot , zrhol(klon) + REAL zchau ,zfroi ,zfice(klon),zneb(klon) + ! + LOGICAL appel1er + SAVE appel1er + !$OMP THREADPRIVATE(appel1er) + ! + !--------------------------------------------------------------- + ! + !AA Variables traceurs: + !AA Provisoire !!! Parametres alpha du lessivage + !AA A priori on a 4 scavenging # possibles + ! + REAL a_tr_sca(4) + save a_tr_sca + !$OMP THREADPRIVATE(a_tr_sca) + ! + ! Variables intermediaires + ! + REAL zalpha_tr + REAL zfrac_lessi + REAL zprec_cond(klon) + !AA + REAL zmair, zcpair, zcpeau + ! Pour la conversion eau-neige + REAL zlh_solid(klon), zm_solid + !IM + !ym INTEGER klevm1 + !--------------------------------------------------------------- + ! + ! Fonctions en ligne: + ! + REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace + REAL zzz + include "YOETHF.h" + include "FCTTRE.h" + fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con + fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc + ! + DATA appel1er /.TRUE./ + !ym + zdelq=0.0 + + print*,'NUAGES4 A. JAM' + IF (appel1er) THEN + ! + PRINT*, 'fisrtilp, ninter:', ninter + PRINT*, 'fisrtilp, evap_prec:', evap_prec + PRINT*, 'fisrtilp, cpartiel:', cpartiel + IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN + PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime + PRINT*, 'Je prefere un sous-intervalle de 6 minutes' + ! CALL abort + ENDIF + appel1er = .FALSE. + ! + !AA initialiation provisoire + a_tr_sca(1) = -0.5 + a_tr_sca(2) = -0.5 + a_tr_sca(3) = -0.5 + a_tr_sca(4) = -0.5 + ! + !AA Initialisation a 1 des coefs des fractions lessivees + ! + !cdir collapse + DO k = 1, klev + DO i = 1, klon + pfrac_nucl(i,k)=1. + pfrac_1nucl(i,k)=1. + pfrac_impa(i,k)=1. + ENDDO + ENDDO + + ENDIF ! test sur appel1er + ! + !MAf Initialisation a 0 de zoliq + ! DO i = 1, klon + ! zoliq(i)=0. + ! ENDDO + ! Determiner les nuages froids par leur temperature + ! nexpo regle la raideur de la transition eau liquide / eau glace. + ! + ztglace = RTT - 15.0 + nexpo = 6 + !cc nexpo = 1 + ! + ! Initialiser les sorties: + ! + !cdir collapse + DO k = 1, klev+1 + DO i = 1, klon + prfl(i,k) = 0.0 + psfl(i,k) = 0.0 + ENDDO + ENDDO + + !cdir collapse + DO k = 1, klev + DO i = 1, klon + d_t(i,k) = 0.0 + d_q(i,k) = 0.0 + d_ql(i,k) = 0.0 + rneb(i,k) = 0.0 + radliq(i,k) = 0.0 + frac_nucl(i,k) = 1. + frac_impa(i,k) = 1. + ENDDO + ENDDO + DO i = 1, klon + rain(i) = 0.0 + snow(i) = 0.0 + zoliq(i)=0. + ! ENDDO + ! + ! Initialiser le flux de precipitation a zero + ! + ! DO i = 1, klon + zrfl(i) = 0.0 + zneb(i) = seuil_neb + ENDDO + ! + ! + !AA Pour plus de securite + + zalpha_tr = 0. + zfrac_lessi = 0. + + !AA---------------------------------------------------------- + ! + ncoreczq=0 + ! Boucle verticale (du haut vers le bas) + ! + !IM : klevm1 + !ym klevm1=klev-1 + DO k = klev, 1, -1 + ! + !AA---------------------------------------------------------- + ! + DO i = 1, klon + zt(i)=t(i,k) + zq(i)=q(i,k) + ENDDO + ! + ! Calculer la varition de temp. de l'air du a la chaleur sensible + ! transporter par la pluie. + ! Il resterait a rajouter cet effet de la chaleur sensible sur les + ! flux de surface, du a la diff. de temp. entre le 1er niveau et la + ! surface. + ! + IF(k.LE.klevm1) THEN + DO i = 1, klon + !IM + zmair=(paprs(i,k)-paprs(i,k+1))/RG + zcpair=RCPD*(1.0+RVTMP2*zq(i)) + zcpeau=RCPD*RVTMP2 + zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau & + + zmair*zcpair*zt(i) ) & + / (zmair*zcpair + zrfl(i)*dtime*zcpeau) + ! C WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1)) + ENDDO + ENDIF + ! + ! + ! Calculer l'evaporation de la precipitation + ! + + + IF (evap_prec) THEN + DO i = 1, klon + IF (zrfl(i) .GT.0.) THEN + IF (thermcep) THEN + zdelta=MAX(0.,SIGN(1.,RTT-zt(i))) + zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k) + zqs(i)=MIN(0.5,zqs(i)) + zcor=1./(1.-RETV*zqs(i)) + zqs(i)=zqs(i)*zcor + ELSE + IF (zt(i) .LT. t_coup) THEN + zqs(i) = qsats(zt(i)) / pplay(i,k) + ELSE + zqs(i) = qsatl(zt(i)) / pplay(i,k) + ENDIF + ENDIF + zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) ) + zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) & + * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG + zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & + * RG*dtime/(paprs(i,k)-paprs(i,k+1)) + zqev = MIN (zqev, zqevt) + zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & + /RG/dtime + + ! pour la glace, on ré-évapore toute la précip dans la + ! couche du dessous + ! la glace venant de la couche du dessus est simplement + ! dans la couche du dessous. + + IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0. + + zq(i) = zq(i) - (zrfln(i)-zrfl(i)) & + * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime + zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & + * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & + * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) + zrfl(i) = zrfln(i) + ENDIF + ENDDO + ENDIF + ! + ! Calculer Qs et L/Cp*dQs/dT: + ! + IF (thermcep) THEN + DO i = 1, klon + zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) + zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta + zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i)) + zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k) + zqs(i) = MIN(0.5,zqs(i)) + zcor = 1./(1.-RETV*zqs(i)) + zqs(i) = zqs(i)*zcor + zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor) + ENDDO + ELSE + DO i = 1, klon + IF (zt(i).LT.t_coup) THEN + zqs(i) = qsats(zt(i))/pplay(i,k) + zdqs(i) = dqsats(zt(i),zqs(i)) + ELSE + zqs(i) = qsatl(zt(i))/pplay(i,k) + zdqs(i) = dqsatl(zt(i),zqs(i)) + ENDIF + ENDDO + ENDIF + ! + ! Determiner la condensation partielle et calculer la quantite + ! de l'eau condensee: + ! + + IF (cpartiel) THEN + + ! print*,'Dans partiel k=',k + ! + ! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau + ! nuageuse a partir des PDF de Sandrine Bony. + ! rneb : fraction nuageuse + ! zqn : eau totale dans le nuage + ! zcond : eau condensee moyenne dans la maille. + ! on prend en compte le réchauffement qui diminue la partie + ! condensee + ! + ! Version avec les raqts + + if (iflag_pdf.eq.0) then + + do i=1,klon + zdelq = min(ratqs(i,k),0.99) * zq(i) + rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq) + zqn(i) = (zq(i)+zdelq+zqs(i))/2.0 + enddo + + else + ! + ! Version avec les nouvelles PDFs. + do i=1,klon + if(zq(i).lt.1.e-15) then + ncoreczq=ncoreczq+1 + zq(i)=1.e-15 + endif + enddo + + if (iflag_cldcon>=5) then + + call cloudth(klon,klev,k,ztv, & + zq,zqta,fraca, & + qcloud,ctot,zpspsk,paprs,ztla,zthl, & + ratqs,zqs,t) + + do i=1,klon + rneb(i,k)=ctot(i,k) + zqn(i)=qcloud(i) + enddo + + endif + + if (iflag_cldcon <= 4) then + lognormale = .true. + elseif (iflag_cldcon == 6) then + ! lognormale en l'absence des thermiques + lognormale = fraca(:,k) < 1e-10 + else + ! Dans le cas iflag_cldcon=5, on prend systématiquement la + ! bi-gaussienne + lognormale = .false. + end if + + do i=1,klon + if (lognormale(i)) then + zpdf_sig(i)=ratqs(i,k)*zq(i) + zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2)) + zpdf_delta(i)=log(zq(i)/zqs(i)) + zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.)) + zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.)) + zpdf_e1(i)=zpdf_a(i)-zpdf_b(i) + zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i)) + zpdf_e1(i)=1.-erf(zpdf_e1(i)) + zpdf_e2(i)=zpdf_a(i)+zpdf_b(i) + zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i)) + zpdf_e2(i)=1.-erf(zpdf_e2(i)) + endif + enddo + + do i=1,klon + if (lognormale(i)) then + if (zpdf_e1(i).lt.1.e-10) then + rneb(i,k)=0. + zqn(i)=zqs(i) + else + rneb(i,k)=0.5*zpdf_e1(i) + zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i) + endif + endif + + enddo + + + endif ! iflag_pdf + + DO i=1,klon + IF (rneb(i,k) .LE. 0.0) THEN + zqn(i) = 0.0 + rneb(i,k) = 0.0 + zcond(i) = 0.0 + rhcl(i,k)=zq(i)/zqs(i) + ELSE IF (rneb(i,k) .GE. 1.0) THEN + zqn(i) = zq(i) + rneb(i,k) = 1.0 + zcond(i) = MAX(0.0,zqn(i)-zqs(i)) + rhcl(i,k)=1.0 + ELSE + zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k) + rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i) + ENDIF + ENDDO + ! do i=1,klon + ! IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0 + ! IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i) + ! rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k))) + !c zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i)) + !c On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par + !c la convection. + !c ATTENTION !!! Il va falloir verifier tout ca. + ! zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k) + !c print*,'ZDQS ',zdqs(i) + !c--Olivier + ! rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i) + ! IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i) + ! IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0 + !c--fin + ! ENDDO + ELSE + DO i = 1, klon + IF (zq(i).GT.zqs(i)) THEN + rneb(i,k) = 1.0 + ELSE + rneb(i,k) = 0.0 + ENDIF + zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i)) + ENDDO + ENDIF + ! + DO i = 1, klon + zq(i) = zq(i) - zcond(i) + ! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD + zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) + ENDDO + ! + ! Partager l'eau condensee en precipitation et eau liquide nuageuse + ! + DO i = 1, klon + IF (rneb(i,k).GT.0.0) THEN + zoliq(i) = zcond(i) + zrho(i) = pplay(i,k) / zt(i) / RD + zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG) + zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace) + zfice(i) = MIN(MAX(zfice(i),0.0),1.0) + zfice(i) = zfice(i)**nexpo + zneb(i) = MAX(rneb(i,k), seuil_neb) + radliq(i,k) = zoliq(i)/REAL(ninter+1) + ENDIF + ENDDO + ! + DO n = 1, ninter + DO i = 1, klon + IF (rneb(i,k).GT.0.0) THEN + zrhol(i) = zrho(i) * zoliq(i) / zneb(i) + + IF (zneb(i).EQ.seuil_neb) THEN + ztot = 0.0 + ELSE + ! quantite d'eau a eliminer: zchau + ! meme chose pour la glace: zfroi + if (ptconv(i,k)) then + zcl =cld_lc_con + zct =1./cld_tau_con + zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & + *fallvc(zrhol(i)) * zfice(i) + else + zcl =cld_lc_lsc + zct =1./cld_tau_lsc + zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & + *fallvs(zrhol(i)) * zfice(i) + endif + zchau = zct *dtime/REAL(ninter) * zoliq(i) & + *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i)) + ztot = zchau + zfroi + ztot = MAX(ztot ,0.0) + ENDIF + ztot = MIN(ztot,zoliq(i)) + zoliq(i) = MAX(zoliq(i)-ztot , 0.0) + radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1) + ENDIF + ENDDO + ENDDO + ! + DO i = 1, klon + IF (rneb(i,k).GT.0.0) THEN + d_ql(i,k) = zoliq(i) + zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) & + * (paprs(i,k)-paprs(i,k+1))/(RG*dtime) + ENDIF + IF (zt(i).LT.RTT) THEN + psfl(i,k)=zrfl(i) + ELSE + prfl(i,k)=zrfl(i) + ENDIF + ENDDO + ! + ! Calculer les tendances de q et de t: + ! + DO i = 1, klon + d_q(i,k) = zq(i) - q(i,k) + d_t(i,k) = zt(i) - t(i,k) + ENDDO + ! + !AA--------------- Calcul du lessivage stratiforme ------------- + + DO i = 1,klon + ! + zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) & + * (paprs(i,k)-paprs(i,k+1))/RG + IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN + !AA lessivage nucleation LMD5 dans la couche elle-meme + if (t(i,k) .GE. ztglace) THEN + zalpha_tr = a_tr_sca(3) + else + zalpha_tr = a_tr_sca(4) + endif + zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) + pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi) + frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi + ! + ! nucleation avec un facteur -1 au lieu de -0.5 + zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i)) + pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi) + ENDIF + ! + ENDDO ! boucle sur i + ! + !AA Lessivage par impaction dans les couches en-dessous + DO kk = k-1, 1, -1 + DO i = 1, klon + IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN + if (t(i,kk) .GE. ztglace) THEN + zalpha_tr = a_tr_sca(1) + else + zalpha_tr = a_tr_sca(2) + endif + zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) + pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi) + frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi + ENDIF + ENDDO + ENDDO + ! + !AA---------------------------------------------------------- + ! FIN DE BOUCLE SUR K + end DO + ! + !AA----------------------------------------------------------- + ! + ! Pluie ou neige au sol selon la temperature de la 1ere couche + ! + DO i = 1, klon + IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN + snow(i) = zrfl(i) + zlh_solid(i) = RLSTT-RLVTT + ELSE + rain(i) = zrfl(i) + zlh_solid(i) = 0. + ENDIF + ENDDO + ! + ! For energy conservation : when snow is present, the solification + ! latent heat is considered. + DO k = 1, klev + DO i = 1, klon + zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k))) + zmair=(paprs(i,k)-paprs(i,k+1))/RG + zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime + d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair) + END DO + END DO + ! + + if (ncoreczq>0) then + print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.' + endif + +END SUBROUTINE fisrtilp