- Timestamp:
- Mar 24, 2010, 1:41:35 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_env.F90
r970 r1330 1 1 SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, & 2 & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk, zqsat,lev_out)2 & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out) 3 3 4 4 !-------------------------------------------------------------- … … 31 31 REAL zu(ngrid,nlay) 32 32 REAL zv(ngrid,nlay) 33 REAL zqsat(ngrid,nlay)33 REAL pqsat(ngrid,nlay) 34 34 35 INTEGER ig,l ,ll35 INTEGER ig,ll 36 36 37 real zcor,zdelta,zcvm5,qlbef 38 real Tbef,qsatbef 39 real dqsat_dT,DT,num,denom 40 REAL RLvCp,DDT0 41 PARAMETER (DDT0=.01) 42 LOGICAL Zsat 37 real dqsat_dT 38 real RLvCp 43 39 44 Zsat=.false. 45 RLvCp = RLVTT/RCPD 40 logical mask(ngrid,nlay) 46 41 47 ! 48 ! Pr Tprec=Tl calcul de qsat 49 ! Si qsat>qT T=Tl, q=qT 50 ! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt) 51 ! On cherche DDT < DDT0 42 43 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 44 ! Initialisations : 45 !------------------ 46 47 print*,'THERMCELL ENV OPTIMISE ' 48 mask(:,:)=.true. 49 RLvCp = RLVTT/RCPD 50 52 51 ! 53 52 ! calcul des caracteristiques de l environnement … … 57 56 zl(ig,ll)=0. 58 57 zh(ig,ll)=pt(ig,ll) 59 zqsat(ig,ll)=0.60 58 EndDO 61 59 EndDO 62 60 ! 63 61 ! 64 !recherche de saturation dans l environnement 65 DO ll=1,nlay 66 ! les points insatures sont definitifs 67 DO ig=1,ngrid 68 Tbef=pt(ig,ll) 69 zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) 70 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll) 71 qsatbef=MIN(0.5,qsatbef) 72 zcor=1./(1.-retv*qsatbef) 73 qsatbef=qsatbef*zcor 74 Zsat = (max(0.,po(ig,ll)-qsatbef) .gt. 1.e-10) 75 if (Zsat) then 76 qlbef=max(0.,po(ig,ll)-qsatbef) 77 ! si sature: ql est surestime, d'ou la sous-relax 78 DT = 0.5*RLvCp*qlbef 79 ! on pourra enchainer 2 ou 3 calculs sans Do while 80 do while (abs(DT).gt.DDT0) 81 ! il faut verifier si c,a conserve quand on repasse en insature ... 82 Tbef=Tbef+DT 83 zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) 84 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll) 85 qsatbef=MIN(0.5,qsatbef) 86 zcor=1./(1.-retv*qsatbef) 87 qsatbef=qsatbef*zcor 88 ! on veut le signe de qlbef 89 qlbef=po(ig,ll)-qsatbef 90 ! dqsat_dT 91 zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) 92 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta 93 zcor=1./(1.-retv*qsatbef) 94 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor) 95 num=-Tbef+pt(ig,ll)+RLvCp*qlbef 96 denom=1.+RLvCp*dqsat_dT 97 if (denom.lt.1.e-10) then 98 print*,'pb denom' 99 endif 100 DT=num/denom 101 enddo 102 ! on ecrit de maniere conservative (sat ou non) 103 zl(ig,ll) = max(0.,qlbef) 104 ! T = Tl +Lv/Cp ql 105 zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll) 106 zo(ig,ll) = po(ig,ll)-zl(ig,ll) 107 endif 108 !on ecrit zqsat 109 zqsat(ig,ll)=qsatbef 110 EndDO 111 EndDO 62 ! Condensation : 63 !--------------- 64 ! Calcul de l'humidite a saturation et de la condensation 65 66 call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat) 67 DO ll=1,nlay 68 DO ig=1,ngrid 69 zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll)) 70 zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll) ! T = Tl + Lv/Cp ql 71 zo(ig,ll) = po(ig,ll)-zl(ig,ll) 72 ENDDO 73 ENDDO 112 74 ! 113 75 ! 114 76 !----------------------------------------------------------------------- 115 ! incrementation eventuelle de tendances precedentes:116 ! ---------------------------------------------------117 77 118 78 if (prt_level.ge.1) print*,'0 OK convect8' 119 79 120 DO 1010l=1,nlay121 DO 1015ig=1,ngrid122 zpspsk(ig,l )=(pplay(ig,l)/100000.)**RKAPPA123 zu(ig,l )=pu(ig,l)124 zv(ig,l )=pv(ig,l)80 DO ll=1,nlay 81 DO ig=1,ngrid 82 zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA 83 zu(ig,ll)=pu(ig,ll) 84 zv(ig,ll)=pv(ig,ll) 125 85 !attention zh est maintenant le profil de T et plus le profil de theta ! 86 ! Quelle horreur ! A eviter. 126 87 ! 127 88 ! T-> Theta 128 ztv(ig,l )=zh(ig,l)/zpspsk(ig,l)89 ztv(ig,ll)=zh(ig,ll)/zpspsk(ig,ll) 129 90 !Theta_v 130 ztv(ig,l)=ztv(ig,l)*(1.+RETV*(zo(ig,l)) & 131 & -zl(ig,l)) 91 ztv(ig,ll)=ztv(ig,ll)*(1.+RETV*(zo(ig,ll))-zl(ig,ll)) 132 92 !Thetal 133 zthl(ig,l )=pt(ig,l)/zpspsk(ig,l)93 zthl(ig,ll)=pt(ig,ll)/zpspsk(ig,ll) 134 94 ! 135 1015 CONTINUE 136 1010 CONTINUE 95 ENDDO 96 ENDDO 137 97 138 98 RETURN
Note: See TracChangeset
for help on using the changeset viewer.