| 1 | SUBROUTINE cloudth(ngrid,klev,ind2, & |
|---|
| 2 | & ztv,po,zqta,fraca, & |
|---|
| 3 | & qcloud,ctot,zpspsk,paprs,ztla,zthl, & |
|---|
| 4 | & ratqs,zqs,t) |
|---|
| 5 | |
|---|
| 6 | |
|---|
| 7 | IMPLICIT NONE |
|---|
| 8 | |
|---|
| 9 | |
|---|
| 10 | !=========================================================================== |
|---|
| 11 | ! Auteur : Arnaud Octavio Jam (LMD/CNRS) |
|---|
| 12 | ! Date : 25 Mai 2010 |
|---|
| 13 | ! Objet : calcule les valeurs de qc et rneb dans les thermiques |
|---|
| 14 | !=========================================================================== |
|---|
| 15 | |
|---|
| 16 | |
|---|
| 17 | #include "YOMCST.h" |
|---|
| 18 | #include "YOETHF.h" |
|---|
| 19 | #include "FCTTRE.h" |
|---|
| 20 | #include "iniprint.h" |
|---|
| 21 | #include "thermcell.h" |
|---|
| 22 | |
|---|
| 23 | INTEGER itap,ind1,ind2 |
|---|
| 24 | INTEGER ngrid,klev,klon,l,ig |
|---|
| 25 | |
|---|
| 26 | REAL ztv(ngrid,klev) |
|---|
| 27 | REAL po(ngrid) |
|---|
| 28 | REAL zqenv(ngrid) |
|---|
| 29 | REAL zqta(ngrid,klev) |
|---|
| 30 | |
|---|
| 31 | REAL fraca(ngrid,klev+1) |
|---|
| 32 | REAL zpspsk(ngrid,klev) |
|---|
| 33 | REAL paprs(ngrid,klev+1) |
|---|
| 34 | REAL ztla(ngrid,klev) |
|---|
| 35 | REAL zthl(ngrid,klev) |
|---|
| 36 | |
|---|
| 37 | REAL zqsatth(ngrid,klev) |
|---|
| 38 | REAL zqsatenv(ngrid,klev) |
|---|
| 39 | |
|---|
| 40 | |
|---|
| 41 | REAL sigma1(ngrid,klev) |
|---|
| 42 | REAL sigma2(ngrid,klev) |
|---|
| 43 | REAL qlth(ngrid,klev) |
|---|
| 44 | REAL qlenv(ngrid,klev) |
|---|
| 45 | REAL qltot(ngrid,klev) |
|---|
| 46 | REAL cth(ngrid,klev) |
|---|
| 47 | REAL cenv(ngrid,klev) |
|---|
| 48 | REAL ctot(ngrid,klev) |
|---|
| 49 | REAL rneb(ngrid,klev) |
|---|
| 50 | REAL t(ngrid,klev) |
|---|
| 51 | REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi |
|---|
| 52 | REAL rdd,cppd,Lv |
|---|
| 53 | REAL alth,alenv,ath,aenv |
|---|
| 54 | REAL sth,senv,sigma1s,sigma2s,xth,xenv |
|---|
| 55 | REAL Tbef,zdelta,qsatbef,zcor |
|---|
| 56 | REAL alpha,qlbef |
|---|
| 57 | REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur |
|---|
| 58 | |
|---|
| 59 | REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid) |
|---|
| 60 | REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) |
|---|
| 61 | REAL zqs(ngrid), qcloud(ngrid) |
|---|
| 62 | REAL erf |
|---|
| 63 | |
|---|
| 64 | |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | |
|---|
| 68 | ! print*,ngrid,klev,ind1,ind2,ztv(ind1,ind2),po(ind1),zqta(ind1,ind2), & |
|---|
| 69 | ! & fraca(ind1,ind2),zpspsk(ind1,ind2),paprs(ind1,ind2),ztla(ind1,ind2),zthl(ind1,ind2), & |
|---|
| 70 | ! & 'verif' |
|---|
| 71 | |
|---|
| 72 | |
|---|
| 73 | ! LOGICAL active(ngrid) |
|---|
| 74 | |
|---|
| 75 | !----------------------------------------------------------------------------------------------------------------- |
|---|
| 76 | ! Initialisation des variables réelles |
|---|
| 77 | !----------------------------------------------------------------------------------------------------------------- |
|---|
| 78 | sigma1(:,:)=0. |
|---|
| 79 | sigma2(:,:)=0. |
|---|
| 80 | qlth(:,:)=0. |
|---|
| 81 | qlenv(:,:)=0. |
|---|
| 82 | qltot(:,:)=0. |
|---|
| 83 | rneb(:,:)=0. |
|---|
| 84 | qcloud(:)=0. |
|---|
| 85 | cth(:,:)=0. |
|---|
| 86 | cenv(:,:)=0. |
|---|
| 87 | ctot(:,:)=0. |
|---|
| 88 | qsatmmussig1=0. |
|---|
| 89 | qsatmmussig2=0. |
|---|
| 90 | rdd=287.04 |
|---|
| 91 | cppd=1005.7 |
|---|
| 92 | pi=3.14159 |
|---|
| 93 | Lv=2.5e6 |
|---|
| 94 | sqrt2pi=sqrt(2.*pi) |
|---|
| 95 | |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | !------------------------------------------------------------------------------------------------------------------ |
|---|
| 99 | ! Calcul de la fraction du thermique et des écart-types des distributions |
|---|
| 100 | !------------------------------------------------------------------------------------------------------------------ |
|---|
| 101 | do ind1=1,ngrid |
|---|
| 102 | |
|---|
| 103 | if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then |
|---|
| 104 | |
|---|
| 105 | zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) |
|---|
| 106 | |
|---|
| 107 | |
|---|
| 108 | ! zqenv(ind1)=po(ind1) |
|---|
| 109 | Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2) |
|---|
| 110 | zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
|---|
| 111 | qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
|---|
| 112 | qsatbef=MIN(0.5,qsatbef) |
|---|
| 113 | zcor=1./(1.-retv*qsatbef) |
|---|
| 114 | qsatbef=qsatbef*zcor |
|---|
| 115 | zqsatenv(ind1,ind2)=qsatbef |
|---|
| 116 | |
|---|
| 117 | |
|---|
| 118 | |
|---|
| 119 | |
|---|
| 120 | alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) |
|---|
| 121 | aenv=1./(1.+(alenv*Lv/cppd)) |
|---|
| 122 | senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) |
|---|
| 123 | |
|---|
| 124 | |
|---|
| 125 | |
|---|
| 126 | |
|---|
| 127 | Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2) |
|---|
| 128 | zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
|---|
| 129 | qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
|---|
| 130 | qsatbef=MIN(0.5,qsatbef) |
|---|
| 131 | zcor=1./(1.-retv*qsatbef) |
|---|
| 132 | qsatbef=qsatbef*zcor |
|---|
| 133 | zqsatth(ind1,ind2)=qsatbef |
|---|
| 134 | |
|---|
| 135 | alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) |
|---|
| 136 | ath=1./(1.+(alth*Lv/cppd)) |
|---|
| 137 | sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) |
|---|
| 138 | |
|---|
| 139 | |
|---|
| 140 | |
|---|
| 141 | !----------------------------------------------------------------------------------------------------------------- |
|---|
| 142 | ! Calcul des écart-types pour s |
|---|
| 143 | !----------------------------------------------------------------------------------------------------------------- |
|---|
| 144 | |
|---|
| 145 | sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) |
|---|
| 146 | sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) |
|---|
| 147 | |
|---|
| 148 | |
|---|
| 149 | !----------------------------------------------------------------------------------------------------------------- |
|---|
| 150 | ! Calcul de l'eau condensée et de la couverture nuageuse |
|---|
| 151 | !----------------------------------------------------------------------------------------------------------------- |
|---|
| 152 | sqrt2pi=sqrt(2.*pi) |
|---|
| 153 | xth=sth/(sqrt(2.)*sigma2s) |
|---|
| 154 | xenv=senv/(sqrt(2.)*sigma1s) |
|---|
| 155 | cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) |
|---|
| 156 | cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
|---|
| 157 | ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) |
|---|
| 158 | ! ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*cenv(ind1,ind2) |
|---|
| 159 | |
|---|
| 160 | |
|---|
| 161 | |
|---|
| 162 | qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2)) |
|---|
| 163 | qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) |
|---|
| 164 | qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) |
|---|
| 165 | ! qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2) |
|---|
| 166 | |
|---|
| 167 | |
|---|
| 168 | ! print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha' |
|---|
| 169 | |
|---|
| 170 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 171 | if (ctot(ind1,ind2).lt.1.e-10) then |
|---|
| 172 | ctot(ind1,ind2)=0. |
|---|
| 173 | qcloud(ind1)=zqsatenv(ind1,ind2) |
|---|
| 174 | |
|---|
| 175 | else |
|---|
| 176 | |
|---|
| 177 | ctot(ind1,ind2)=ctot(ind1,ind2) |
|---|
| 178 | qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1) |
|---|
| 179 | |
|---|
| 180 | endif |
|---|
| 181 | |
|---|
| 182 | |
|---|
| 183 | ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' |
|---|
| 184 | |
|---|
| 185 | |
|---|
| 186 | else ! gaussienne environnement seule |
|---|
| 187 | |
|---|
| 188 | zqenv(ind1)=po(ind1) |
|---|
| 189 | Tbef=t(ind1,ind2) |
|---|
| 190 | zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) |
|---|
| 191 | qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2) |
|---|
| 192 | qsatbef=MIN(0.5,qsatbef) |
|---|
| 193 | zcor=1./(1.-retv*qsatbef) |
|---|
| 194 | qsatbef=qsatbef*zcor |
|---|
| 195 | zqsatenv(ind1,ind2)=qsatbef |
|---|
| 196 | |
|---|
| 197 | |
|---|
| 198 | ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) |
|---|
| 199 | zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd) |
|---|
| 200 | alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) |
|---|
| 201 | aenv=1./(1.+(alenv*Lv/cppd)) |
|---|
| 202 | senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) |
|---|
| 203 | |
|---|
| 204 | |
|---|
| 205 | sigma1s=ratqs(ind1,ind2)*zqenv(ind1) |
|---|
| 206 | |
|---|
| 207 | sqrt2pi=sqrt(2.*pi) |
|---|
| 208 | xenv=senv/(sqrt(2.)*sigma1s) |
|---|
| 209 | ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) |
|---|
| 210 | qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2)) |
|---|
| 211 | |
|---|
| 212 | if (ctot(ind1,ind2).lt.1.e-3) then |
|---|
| 213 | ctot(ind1,ind2)=0. |
|---|
| 214 | qcloud(ind1)=zqsatenv(ind1,ind2) |
|---|
| 215 | |
|---|
| 216 | else |
|---|
| 217 | |
|---|
| 218 | ctot(ind1,ind2)=ctot(ind1,ind2) |
|---|
| 219 | qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2) |
|---|
| 220 | |
|---|
| 221 | endif |
|---|
| 222 | |
|---|
| 223 | |
|---|
| 224 | |
|---|
| 225 | |
|---|
| 226 | |
|---|
| 227 | |
|---|
| 228 | endif |
|---|
| 229 | enddo |
|---|
| 230 | |
|---|
| 231 | return |
|---|
| 232 | end |
|---|
| 233 | |
|---|
| 234 | |
|---|
| 235 | |
|---|
| 236 | |
|---|
| 237 | |
|---|
| 238 | |
|---|
| 239 | |
|---|
| 240 | |
|---|
| 241 | |
|---|
| 242 | |
|---|
| 243 | |
|---|
