Changeset 974
- Timestamp:
- Jun 19, 2008, 12:26:15 PM (17 years ago)
- Location:
- LMDZ4/trunk/libf/phylmd
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/trunk/libf/phylmd/calwake.F
r953 r974 30 30 IMPLICIT none 31 31 c====================================================================== 32 #include "dimensions.h" 33 !#include "dimphy.h" 34 #include "YOMCST.h" 35 36 c Arguments 37 c---------- 38 39 INTEGER i,l,ktopw(klon) 40 REAL dtime 41 42 REAL paprs(klon,klev+1),pplay(klon,klev) 43 REAL t(klon,klev), q(klon,klev), omgb(klon,klev) 44 REAL dt_dwn(klon,klev), dq_dwn(klon,klev),M_dwn(klon,klev) 45 REAL M_up(klon,klev) 46 REAL dt_a(klon,klev), dq_a(klon,klev) 47 REAL wdt_PBL(klon,klev), wdq_PBL(klon,klev) 48 REAL udt_PBL(klon,klev), udq_PBL(klon,klev) 49 REAL wake_deltat(klon,klev),wake_deltaq(klon,klev) 50 REAL dt_wake(klon,klev),dq_wake(klon,klev) 51 REAL wake_d_deltat_gw(klon,klev) 52 REAL wake_h(klon),wake_s(klon) 53 REAL wake_dth(klon,klev) 54 REAL wake_pe(klon),wake_fip(klon),wake_gfl(klon) 55 REAL undi_t(klon,klev),undi_q(klon,klev) 56 REAL wake_omgbdth(klon,klev),wake_dp_omgb(klon,klev) 57 REAL wake_dtKE(klon,klev),wake_dqKE(klon,klev) 58 REAL wake_dtPBL(klon,klev),wake_dqPBL(klon,klev) 59 REAL wake_omg(klon,klev+1),wake_dp_deltomg(klon,klev) 60 REAL wake_spread(klon,klev),wake_Cstar(klon) 61 REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev) 62 REAL d_deltatw(klon,klev), d_deltaqw(klon,klev) 63 INTEGER wake_k(klon) 64 REAL sigd(klon) 65 REAL wake_dens(klon) 66 67 C Variable internes 68 C ----------------- 69 70 REAL aire 71 REAL p(klon,klev),ph(klon,klev+1),pi(klon,klev) 72 REAL te(klon,klev),qe(klon,klev),omgbe(klon,klev) 73 REAL dtdwn(klon,klev),dqdwn(klon,klev) 74 REAL dta(klon,klev),dqa(klon,klev) 75 REAL wdtPBL(klon,klev),wdqPBL(klon,klev) 76 REAL udtPBL(klon,klev),udqPBL(klon,klev) 77 REAL amdwn(klon,klev),amup(klon,klev) 78 REAL dtw(klon,klev),dqw(klon,klev),dth(klon,klev) 79 REAL d_deltat_gw(klon,klev) 80 REAL dtls(klon,klev),dqls(klon,klev) 81 REAL tu(klon,klev),qu(klon,klev) 82 REAL hw(klon),sigmaw(klon),wape(klon),fip(klon),gfl(klon) 83 REAL omgbdth(klon,klev),dp_omgb(klon,klev) 84 REAL dtKE(klon,klev),dqKE(klon,klev) 85 REAL dtPBL(klon,klev),dqPBL(klon,klev) 86 REAL omg(klon,klev+1),dp_deltomg(klon,klev),spread(klon,klev) 87 REAL Cstar(klon) 88 REAL sigd0(klon),wdens(klon) 89 90 REAL RDCP 91 92 c print *, '-> calwake, wake_s ', wake_s(1) 93 94 RDCP=1./3.5 95 96 c----------------------------------------------------------- 97 cIM 290108 DO 999 i=1,klon ! a vectoriser 98 c---------------------------------------------------------- 99 100 101 DO l=1,klev 102 DO i=1,klon 103 p(i,l)= pplay(i,l) 104 ph(i,l)= paprs(i,l) 105 pi(i,l) = (pplay(i,l)/100000.)**RDCP 106 107 te(i,l) = t(i,l) 108 qe(i,l) = q(i,l) 109 omgbe(i,l) = omgb(i,l) 110 111 dtdwn(i,l)= dt_dwn(i,l) 112 dqdwn(i,l)= dq_dwn(i,l) 113 dta(i,l)= dt_a(i,l) 114 dqa(i,l)= dq_a(i,l) 115 wdtPBL(i,l)= wdt_PBL(i,l) 116 wdqPBL(i,l)= wdq_PBL(i,l) 117 udtPBL(i,l)= udt_PBL(i,l) 118 udqPBL(i,l)= udq_PBL(i,l) 119 ENDDO 120 ENDDO 121 122 DO i=1,klon 123 sigd0(i)=sigd(i) 124 ENDDO 125 c print*, 'sigd0,sigd', sigd0, sigd(i) 126 DO i=1,klon 127 ph(i,klev+1)=0. 128 ENDDO 129 130 DO i=1,klon 131 ktopw(i) = wake_k(i) 132 ENDDO 133 134 DO l=1,klev 135 DO i=1,klon 136 dtw(i,l) = wake_deltat(i,l) 137 dqw(i,l) = wake_deltaq(i,l) 138 ENDDO 139 ENDDO 140 141 DO l=1,klev 142 DO i=1,klon 143 dtls(i,l)=dt_wake(i,l) 144 dqls(i,l)=dq_wake(i,l) 145 ENDDO 146 ENDDO 147 148 DO i=1,klon 149 hw(i) = wake_h(i) 150 sigmaw(i)= wake_s(i) 151 ENDDO 152 153 cfkc les flux de masses sont evalues aux niveaux et valent 0 a la surface 154 cfkc on veut le flux de masse au milieu des couches 155 156 DO l=1,klev-1 157 DO i=1,klon 158 amdwn(i,l)= 0.5*(M_dwn(i,l)+M_dwn(i,l+1)) 159 amdwn(i,l)= (M_dwn(i,l+1)) 160 ENDDO 161 ENDDO 162 163 c au sommet le flux de masse est nul 164 165 DO i=1,klon 166 amdwn(i,klev)=0.5*M_dwn(i,klev) 167 ENDDO 168 c 169 DO l = 1,klev 170 DO i=1,klon 171 amup(i,l)=M_up(i,l) 172 ENDDO 173 ENDDO 174 175 call WAKE(p,ph,pi,dtime,sigd0 176 $ ,te,qe,omgbe 177 $ ,dtdwn,dqdwn,amdwn,amup,dta,dqa 178 $ ,wdtPBL,wdqPBL,udtPBL,udqPBL 179 $ ,dtw,dqw,dth,hw,sigmaw,wape,fip,gfl 180 $ ,dtls,dqls,ktopw 181 $ ,omgbdth,dp_omgb,wdens 182 $ ,tu,qu 183 $ ,dtKE,dqKE 184 $ ,dtPBL,dqPBL 185 $ ,omg,dp_deltomg,spread 186 $ ,Cstar,d_deltat_gw 187 $ ,d_deltatw,d_deltaqw) 188 189 DO i=1,klon 190 IF (ktopw(i) .GT. 0) THEN 191 DO l=1,klev 192 wake_deltat(i,l)= dtw(i,l) 193 wake_deltaq(i,l)= dqw(i,l) 194 wake_d_deltat_gw(i,l)= d_deltat_gw(i,l) 195 wake_omgbdth(i,l) = omgbdth(i,l) 196 wake_dp_omgb(i,l) = dp_omgb(i,l) 197 wake_dtKE(i,l) = dtKE(i,l) 198 wake_dqKE(i,l) = dqKE(i,l) 199 wake_dtPBL(i,l) = dtPBL(i,l) 200 wake_dqPBL(i,l) = dqPBL(i,l) 201 wake_omg(i,l) = omg(i,l) 202 wake_dp_deltomg(i,l) = dp_deltomg(i,l) 203 wake_spread(i,l) = spread(i,l) 204 wake_dth(i,l) = dth(i,l) 205 dt_wake(i,l) = dtls(i,l) 206 dq_wake(i,l) = dqls(i,l) 207 undi_t(i,l) = tu(i,l) 208 undi_q(i,l) = qu(i,l) 209 wake_ddeltat(i,l) = d_deltatw(i,l) 210 wake_ddeltaq(i,l) = d_deltaqw(i,l) 211 ENDDO 212 ELSE 213 DO l = 1,klev 214 wake_deltat(i,l)= 0. 215 wake_deltaq(i,l)= 0. 216 wake_d_deltat_gw(i,l)= 0. 217 wake_omgbdth(i,l) = 0. 218 wake_dp_omgb(i,l) = 0. 219 wake_dtKE(i,l) = 0. 220 wake_dqKE(i,l) = 0. 221 wake_omg(i,l) = 0. 222 wake_dp_deltomg(i,l) = 0. 223 wake_spread(i,l) = 0. 224 wake_dth(i,l)=0. 225 dt_wake(i,l)=0. 226 dq_wake(i,l)=0. 227 undi_t(i,l)=te(i,l) 228 undi_q(i,l)=qe(i,l) 229 ENDDO 230 ENDIF 231 232 wake_h(i)= hw(i) 233 wake_s(i)= sigmaw(i) 234 wake_pe(i)= wape(i) 235 wake_fip(i)= fip(i) 236 wake_gfl(i) = gfl(i) 237 wake_k(i) =ktopw(i) 238 wake_Cstar(i) = Cstar(i) 239 wake_dens(i) = wdens(i) 240 c 241 cIM 290108 999 CONTINUE 242 c 243 ENDDO 244 RETURN 245 END 246 SUBROUTINE CALWAKE_scal(paprs,pplay,dtime 247 : ,t,q,omgb 248 : ,dt_dwn,dq_dwn,M_dwn,M_up 249 : ,dt_a,dq_a,sigd 250 : ,wdt_PBL,wdq_PBL 251 : ,udt_PBL,udq_PBL 252 o ,wake_deltat,wake_deltaq,wake_dth 253 o ,wake_h,wake_s,wake_dens 254 o ,wake_pe,wake_fip,wake_gfl 255 o ,dt_wake,dq_wake 256 o ,wake_k 257 o ,undi_t,undi_q 258 o ,wake_omgbdth,wake_dp_omgb 259 o ,wake_dtKE,wake_dqKE 260 o ,wake_dtPBL,wake_dqPBL 261 o ,wake_omg,wake_dp_deltomg 262 o ,wake_spread,wake_Cstar,wake_d_deltat_gw 263 o ,wake_ddeltat,wake_ddeltaq) 264 *************************************************************** 265 * * 266 * CALWAKE * 267 * interface avec le schema de calcul de la poche * 268 * froide * 269 * * 270 * written by : CHERUY Frederique, 13/03/2000, 10.31.05 * 271 * modified by : ROEHRIG Romain, 01/30/2007 * 272 *************************************************************** 273 * 274 USE dimphy 275 IMPLICIT none 276 c====================================================================== 32 277 33 278 #include "dimensions.h" … … 151 396 ENDDO 152 397 153 call WAKE (p,ph,pi,dtime,sigd0398 call WAKE_scal(p,ph,pi,dtime,sigd0 154 399 $ ,te,qe,omgbe 155 400 $ ,dtdwn,dqdwn,amdwn,amup,dta,dqa -
LMDZ4/trunk/libf/phylmd/wake.F
r953 r974 22 22 *************************************************************** 23 23 c 24 USEdimphy24 use dimphy 25 25 IMPLICIT none 26 26 c============================================================================ … … 112 112 113 113 #include "dimensions.h" 114 #include "YOMCST.h" 115 #include "cvthermo.h" 116 #include "iniprint.h" 117 118 c Arguments en entree 119 c-------------------- 120 121 REAL, dimension(klon,klev) :: p, ppi 122 REAL, dimension(klon,klev+1) :: ph, omgb 123 REAL dtime 124 REAL, dimension(klon,klev) :: te0,qe0 125 REAL, dimension(klon,klev) :: dtdwn, dqdwn 126 REAL, dimension(klon,klev) :: wdtPBL,wdqPBL 127 REAL, dimension(klon,klev) :: udtPBL,udqPBL 128 REAL, dimension(klon,klev) :: amdwn, amup 129 REAL, dimension(klon,klev) :: dta, dqa 130 REAL, dimension(klon) :: sigd_con 131 132 c Sorties 133 c-------- 134 135 REAL, dimension(klon,klev) :: deltatw, deltaqw, dth 136 REAL, dimension(klon,klev) :: tu, qu 137 REAL, dimension(klon,klev) :: dtls, dqls 138 REAL, dimension(klon,klev) :: dtKE, dqKE 139 REAL, dimension(klon,klev) :: dtPBL, dqPBL 140 REAL, dimension(klon,klev) :: spread 141 REAL, dimension(klon,klev) :: d_deltatgw 142 REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2 143 REAL, dimension(klon,klev+1) :: omgbdth, omg 144 REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg 145 REAL, dimension(klon,klev) :: d_deltat_gw 146 REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar 147 INTEGER, dimension(klon) :: ktopw 148 149 c Variables internes 150 c------------------- 151 152 c Variables à fixer 153 REAL ALON 154 REAL coefgw 155 REAL :: wdens0, wdens 156 REAL stark 157 REAL alpk 158 REAL delta_t_min 159 REAL Pupper 160 INTEGER nsub 161 REAL dtimesub 162 REAL sigmad, hwmin 163 cIM 080208 164 LOGICAL, dimension(klon) :: gwake 165 166 c Variables de sauvegarde 167 REAL, dimension(klon,klev) :: deltatw0 168 REAL, dimension(klon,klev) :: deltaqw0 169 REAL, dimension(klon,klev) :: te, qe 170 REAL, dimension(klon) :: sigmaw0, sigmaw1 171 172 c Variables pour les GW 173 REAL, DIMENSION(klon) :: LL 174 REAL, dimension(klon,klev) :: N2 175 REAL, dimension(klon,klev) :: Cgw 176 REAL, dimension(klon,klev) :: Tgw 177 178 c Variables liées au calcul de hw 179 REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new 180 REAL, DIMENSION(klon) :: sum_dth 181 REAL, DIMENSION(klon) :: dthmin 182 REAL, DIMENSION(klon) :: z, dz, hw0 183 INTEGER, DIMENSION(klon) :: ktop, kupper 184 185 c Autres variables internes 186 INTEGER isubstep, k, i 187 188 REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu 189 REAL, DIMENSION(klon) :: sum_dq, sum_rho 190 REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn 191 REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu 192 REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho 193 REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn 194 195 REAL, DIMENSION(klon,klev) :: rho, rhow 196 REAL, DIMENSION(klon,klev+1) :: rhoh 197 REAL, DIMENSION(klon,klev) :: rhow_moyen 198 REAL, DIMENSION(klon,klev) :: zh 199 REAL, DIMENSION(klon,klev+1) :: zhh 200 REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2 201 202 REAL, DIMENSION(klon,klev) :: the, thu 203 204 REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw 205 206 REAL, DIMENSION(klon,klev+1) :: omgbw 207 REAL, DIMENSION(klon) :: omgtop 208 REAL, DIMENSION(klon,klev) :: dp_omgbw 209 REAL, DIMENSION(klon) :: ztop, dztop 210 REAL, DIMENSION(klon,klev) :: alpha_up 211 212 REAL, dimension(klon) :: RRe1, RRe2 213 REAL :: RRd1, RRd2 214 REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2 215 REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth 216 REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq 217 REAL, DIMENSION(klon,klev) :: omgbdq 218 219 REAL, dimension(klon) :: ff, gg 220 REAL, dimension(klon) :: wape2, Cstar2, heff 221 222 REAL, DIMENSION(klon,klev) :: Crep 223 REAL Crep_upper, Crep_sol 224 225 C------------------------------------------------------------------------- 226 c Initialisations 227 c------------------------------------------------------------------------- 228 229 c print*, 'wake initialisations' 230 231 c Essais d'initialisation avec sigmaw = 0.02 et hw = 10. 232 c------------------------------------------------------------------------- 233 234 DATA sigmad, hwmin /.02,10./ 235 236 C Longueur de maille (en m) 237 c------------------------------------------------------------------------- 238 239 c ALON = 3.e5 240 ALON = 1.e6 241 242 243 C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) 244 c 245 c coefgw : Coefficient pour les ondes de gravité 246 c stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) 247 c wdens : Densité de poche froide par maille 248 c------------------------------------------------------------------------- 249 250 coefgw=10 251 c coefgw=1 252 c wdens0 = 1.0/(alon**2) 253 wdens = 1.0/(alon**2) 254 stark = 0.50 255 cCRtest 256 alpk=0.1 257 c alpk = 1.0 258 c alpk = 0.5 259 c alpk = 0.05 260 Crep_upper=0.9 261 Crep_sol=1.0 262 263 264 C Minimum value for |T_wake - T_undist|. Used for wake top definition 265 c------------------------------------------------------------------------- 266 267 delta_t_min = 0.2 268 269 C 1. - Save initial values and initialize tendencies 270 C -------------------------------------------------- 271 272 DO k=1,klev 273 DO i=1, klon 274 deltatw0(i,k) = deltatw(i,k) 275 deltaqw0(i,k)= deltaqw(i,k) 276 te(i,k) = te0(i,k) 277 qe(i,k) = qe0(i,k) 278 dtls(i,k) = 0. 279 dqls(i,k) = 0. 280 d_deltat_gw(i,k)=0. 281 !IM 060508 beg 282 d_deltatw2(i,k)=0. 283 d_deltaqw2(i,k)=0. 284 !IM 060508 end 285 ENDDO 286 ENDDO 287 c sigmaw1=sigmaw 288 c IF (sigd_con.GT.sigmaw1) THEN 289 c print*, 'sigmaw,sigd_con', sigmaw, sigd_con 290 c ENDIF 291 DO i=1, klon 292 cc sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) 293 sigmaw(i) = amax1(sigmaw(i),sigmad) 294 sigmaw(i) = amin1(sigmaw(i),0.99) 295 sigmaw0(i) = sigmaw(i) 296 ENDDO 297 C 298 C 299 C 2. - Prognostic part 300 C -------------------- 301 C 302 C 303 C 2.1 - Undisturbed area and Wake integrals 304 C --------------------------------------------------------- 305 306 DO i=1, klon 307 z(i) = 0. 308 ktop(i)=0 309 kupper(i) = 0 310 sum_thu(i) = 0. 311 sum_tu(i) = 0. 312 sum_qu(i) = 0. 313 sum_thvu(i) = 0. 314 sum_dth(i) = 0. 315 sum_dq(i) = 0. 316 sum_rho(i) = 0. 317 sum_dtdwn(i) = 0. 318 sum_dqdwn(i) = 0. 319 320 av_thu(i) = 0. 321 av_tu(i) =0. 322 av_qu(i) =0. 323 av_thvu(i) = 0. 324 av_dth(i) = 0. 325 av_dq(i) = 0. 326 av_rho(i) =0. 327 av_dtdwn(i) =0. 328 av_dqdwn(i) = 0. 329 ENDDO 330 c 331 c Distance between wakes 332 DO i = 1,klon 333 LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens) 334 ENDDO 335 C Potential temperatures and humidity 336 c---------------------------------------------------------- 337 DO k =1,klev 338 DO i=1, klon 339 rho(i,k) = p(i,k)/(rd*te(i,k)) 340 IF(k .eq. 1) THEN 341 rhoh(i,k) = ph(i,k)/(rd*te(i,k)) 342 zhh(i,k)=0 343 ELSE 344 rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1))) 345 zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1) 346 ENDIF 347 the(i,k) = te(i,k)/ppi(i,k) 348 thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k) 349 tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i) 350 qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i) 351 rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k))) 352 dth(i,k) = deltatw(i,k)/ppi(i,k) 353 ENDDO 354 ENDDO 355 356 DO k = 1, klev-1 357 DO i=1, klon 358 IF(k.eq.1) THEN 359 N2(i,k)=0 360 ELSE 361 N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)- 362 $ the(i,k-1))/(p(i,k+1)-p(i,k-1))) 363 ENDIF 364 ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2 365 366 Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k) 367 Tgw(i,k)=coefgw*Cgw(i,k)/LL(i) 368 ENDDO 369 ENDDO 370 371 DO i=1, klon 372 N2(i,klev)=0 373 ZH(i,klev)=0 374 Cgw(i,klev)=0 375 Tgw(i,klev)=0 376 ENDDO 377 378 c Calcul de la masse volumique moyenne de la colonne (bdlmd) 379 c----------------------------------------------------------------- 380 381 DO k=1,klev 382 DO i=1, klon 383 epaisseur1(i,k)=0. 384 epaisseur2(i,k)=0. 385 ENDDO 386 ENDDO 387 388 DO i=1, klon 389 epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1. 390 epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1. 391 rhow_moyen(i,1) = rhow(i,1) 392 ENDDO 393 394 DO k = 2, klev 395 DO i=1, klon 396 epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1. 397 epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k) 398 rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+ 399 $ rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k) 400 ENDDO 401 ENDDO 402 403 C 404 C Choose an integration bound well above wake top 405 c----------------------------------------------------------------- 406 c 407 C Pupper = 50000. ! melting level 408 Pupper = 60000. 409 c Pupper = 80000. ! essais pour case_e 410 C 411 C Determine Wake top pressure (Ptop) from buoyancy integral 412 C -------------------------------------------------------- 413 c 414 c-1/ Pressure of the level where dth becomes less than delta_t_min. 415 416 DO i=1,klon 417 ptop_provis(i)=ph(i,1) 418 ENDDO 419 DO k= 2,klev 420 DO i=1,klon 421 c 422 cIM v3JYG; ptop_provis(i).LT. ph(i,1) 423 c 424 IF (dth(i,k) .GT. -delta_t_min .and. 425 $ dth(i,k-1).LT. -delta_t_min .and. 426 $ ptop_provis(i).EQ. ph(i,1)) THEN 427 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) 428 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) / 429 $ (dth(i,k) - dth(i,k-1)) 430 ENDIF 431 ENDDO 432 ENDDO 433 434 c-2/ dth integral 435 436 DO i=1,klon 437 sum_dth(i) = 0. 438 dthmin(i) = -delta_t_min 439 z(i) = 0. 440 ENDDO 441 442 DO k = 1,klev 443 DO i=1,klon 444 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg) 445 IF (dz(i) .gt. 0) THEN 446 z(i) = z(i)+dz(i) 447 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 448 dthmin(i) = amin1(dthmin(i),dth(i,k)) 449 ENDIF 450 ENDDO 451 ENDDO 452 453 c-3/ height of triangle with area= sum_dth and base = dthmin 454 455 DO i=1,klon 456 hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5) 457 hw0(i) = amax1(hwmin,hw0(i)) 458 ENDDO 459 460 c-4/ now, get Ptop 461 462 DO i=1,klon 463 z(i) = 0. 464 ptop(i) = ph(i,1) 465 ENDDO 466 467 DO k = 1,klev 468 DO i=1,klon 469 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i)) 470 IF (dz(i) .gt. 0) THEN 471 z(i) = z(i)+dz(i) 472 ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i) 473 ENDIF 474 ENDDO 475 ENDDO 476 477 478 C-5/ Determination de ktop et kupper 479 480 DO k=klev,1,-1 481 DO i=1,klon 482 IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k 483 IF (ph(i,k+1) .lt. pupper) kupper(i)=k 484 ENDDO 485 ENDDO 486 487 c-6/ Correct ktop and ptop 488 489 DO i = 1,klon 490 ptop_new(i)=ptop(i) 491 ENDDO 492 DO k= klev,2,-1 493 DO i=1,klon 494 IF (k .LE. ktop(i) .and. 495 $ ptop_new(i) .EQ. ptop(i) .and. 496 $ dth(i,k) .GT. -delta_t_min .and. 497 $ dth(i,k-1).LT. -delta_t_min) THEN 498 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) 499 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) / 500 $ (dth(i,k) - dth(i,k-1)) 501 ENDIF 502 ENDDO 503 ENDDO 504 505 DO i=1,klon 506 ptop(i) = ptop_new(i) 507 ENDDO 508 509 DO k=klev,1,-1 510 DO i=1,klon 511 IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k 512 ENDDO 513 ENDDO 514 c 515 c-5/ Set deltatw & deltaqw to 0 above kupper 516 c 517 DO k = 1,klev 518 DO i=1,klon 519 IF (k.GE. kupper(i)) THEN 520 deltatw(i,k) = 0. 521 deltaqw(i,k) = 0. 522 ENDIF 523 ENDDO 524 ENDDO 525 c 526 C 527 C Vertical gradient of LS omega 528 C 529 DO k = 1,klev 530 DO i=1,klon 531 IF (k.LE. kupper(i)) THEN 532 dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k)) 533 ENDIF 534 ENDDO 535 ENDDO 536 C 537 C Integrals (and wake top level number) 538 C -------------------------------------- 539 C 540 C Initialize sum_thvu to 1st level virt. pot. temp. 541 542 DO i=1,klon 543 z(i) = 1. 544 dz(i) = 1. 545 sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) 546 sum_dth(i) = 0. 547 ENDDO 548 549 DO k = 1,klev 550 DO i=1,klon 551 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) 552 IF (dz(i) .GT. 0) THEN 553 z(i) = z(i)+dz(i) 554 sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i) 555 sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i) 556 sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i) 557 sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i) 558 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 559 sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i) 560 sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i) 561 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i) 562 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i) 563 ENDIF 564 ENDDO 565 ENDDO 566 c 567 DO i=1,klon 568 hw0(i) = z(i) 569 ENDDO 570 c 571 C 572 C 2.1 - WAPE and mean forcing computation 573 C --------------------------------------- 574 C 575 C --------------------------------------- 576 C 577 C Means 578 579 DO i=1,klon 580 av_thu(i) = sum_thu(i)/hw0(i) 581 av_tu(i) = sum_tu(i)/hw0(i) 582 av_qu(i) = sum_qu(i)/hw0(i) 583 av_thvu(i) = sum_thvu(i)/hw0(i) 584 c av_thve = sum_thve/hw0 585 av_dth(i) = sum_dth(i)/hw0(i) 586 av_dq(i) = sum_dq(i)/hw0(i) 587 av_rho(i) = sum_rho(i)/hw0(i) 588 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 589 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 590 591 wape(i) = - rg*hw0(i)*(av_dth(i) 592 $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)* 593 $ av_dq(i) ))/av_thvu(i) 594 ENDDO 595 C 596 C 2.2 Prognostic variable update 597 C ------------------------------ 598 C 599 C Filter out bad wakes 600 601 DO k = 1,klev 602 DO i=1,klon 603 IF ( wape(i) .LT. 0.) THEN 604 deltatw(i,k) = 0. 605 deltaqw(i,k) = 0. 606 dth(i,k) = 0. 607 ENDIF 608 ENDDO 609 ENDDO 610 c 611 DO i=1,klon 612 IF ( wape(i) .LT. 0.) THEN 613 wape(i) = 0. 614 Cstar(i) = 0. 615 hw(i) = hwmin 616 sigmaw(i) = amax1(sigmad,sigd_con(i)) 617 fip(i) = 0. 618 gwake(i) = .FALSE. 619 ELSE 620 Cstar(i) = stark*sqrt(2.*wape(i)) 621 gwake(i) = .TRUE. 622 ENDIF 623 ENDDO 624 c 625 C 626 CC ----------------------------------------------------------------- 627 C Sub-time-stepping 628 C ----------------- 629 C 630 nsub=10 631 dtimesub=dtime/nsub 632 c 633 c------------------------------------------------------------ 634 DO isubstep = 1,nsub 635 c------------------------------------------------------------ 636 DO i=1,klon 637 gfl(i) = 2.*sqrt(3.14*wdens*sigmaw(i)) 638 ENDDO 639 DO i=1,klon 640 sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub 641 sigmaw(i) =amin1(sigmaw(i),0.99) !!!!!!!! 642 c wdens = wdens0/(10.*sigmaw) 643 c sigmaw =max(sigmaw,sigd_con) 644 c sigmaw =max(sigmaw,sigmad) 645 ENDDO 646 C 647 C 648 c calcul de la difference de vitesse verticale poche - zone non perturbee 649 cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg 650 cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit 651 cIM 060208 au niveau k=1..? 652 dp_deltomg(1:klon,1:klev)=0. 653 DO k= 1,klev+1 654 DO i = 1,klon 655 omg(i,k)=0. 656 ENDDO 657 ENDDO 658 c 659 DO i=1,klon 660 z(i)= 0. 661 omg(i,1) = 0. 662 dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i))) 663 ENDDO 664 c 665 DO k= 2,klev 666 DO i = 1,klon 667 IF( k .LE. ktop(i)) THEN 668 dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg) 669 z(i) = z(i)+dz(i) 670 dp_deltomg(i,k)= dp_deltomg(i,1) 671 omg(i,k)= dp_deltomg(i,1)*z(i) 672 ENDIF 673 ENDDO 674 ENDDO 675 c 676 DO i = 1,klon 677 dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg) 678 ztop(i) = z(i)+dztop(i) 679 omgtop(i)=dp_deltomg(i,1)*ztop(i) 680 ENDDO 681 c 682 c ----------------- 683 c From m/s to Pa/s 684 c ----------------- 685 c 686 DO i=1,klon 687 omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i) 688 dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1)) 689 ENDDO 690 c 691 DO k= 1,klev 692 DO i = 1,klon 693 IF( k .LE. ktop(i)) THEN 694 omg(i,k) = - rho(i,k)*rg*omg(i,k) 695 dp_deltomg(i,k) = dp_deltomg(i,1) 696 ENDIF 697 ENDDO 698 ENDDO 699 c 700 c raccordement lineaire de omg de ptop a pupper 701 702 DO i=1,klon 703 IF (kupper(i) .GT. ktop(i)) THEN 704 omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i) 705 $ + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i)) 706 dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ 707 $ (ptop(i)-pupper) 708 ENDIF 709 ENDDO 710 c 711 DO k= 1,klev 712 DO i = 1,klon 713 IF( k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN 714 dp_deltomg(i,k) = dp_deltomg(i,kupper(i)) 715 omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i)) 716 ENDIF 717 ENDDO 718 ENDDO 719 c 720 c-- Compute wake average vertical velocity omgbw 721 c 722 c 723 DO k = 1,klev+1 724 DO i=1,klon 725 omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k) 726 ENDDO 727 ENDDO 728 c-- and its vertical gradient dp_omgbw 729 c 730 DO k = 1,klev 731 DO i=1,klon 732 dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k)) 733 ENDDO 734 ENDDO 735 C 736 c-- Upstream coefficients for omgb velocity 737 c-- (alpha_up(k) is the coefficient of the value at level k) 738 c-- (1-alpha_up(k) is the coefficient of the value at level k-1) 739 DO k = 1,klev 740 DO i=1,klon 741 alpha_up(i,k) = 0. 742 IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1. 743 ENDDO 744 ENDDO 745 746 c Matrix expressing [The,deltatw] from [Th1,Th2] 747 748 DO i=1,klon 749 RRe1(i) = 1.-sigmaw(i) 750 RRe2(i) = sigmaw(i) 751 ENDDO 752 RRd1 = -1. 753 RRd2 = 1. 754 c 755 c-- Get [Th1,Th2], dth and [q1,q2] 756 c 757 DO k= 1,klev 758 DO i = 1,klon 759 IF(k .LE. kupper(i)+1) THEN 760 dth(i,k) = deltatw(i,k)/ppi(i,k) 761 Th1(i,k) = the(i,k) - sigmaw(i) *dth(i,k) ! undisturbed area 762 Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k) ! wake 763 q1(i,k) = qe(i,k) - sigmaw(i) *deltaqw(i,k) ! undisturbed area 764 q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake 765 ENDIF 766 ENDDO 767 ENDDO 768 769 DO i=1,klon 770 D_Th1(i,1) = 0. 771 D_Th2(i,1) = 0. 772 D_dth(i,1) = 0. 773 D_q1(i,1) = 0. 774 D_q2(i,1) = 0. 775 D_dq(i,1) = 0. 776 ENDDO 777 778 DO k= 2,klev 779 DO i = 1,klon 780 IF(k .LE. kupper(i)+1) THEN 781 D_Th1(i,k) = Th1(i,k-1)-Th1(i,k) 782 D_Th2(i,k) = Th2(i,k-1)-Th2(i,k) 783 D_dth(i,k) = dth(i,k-1)-dth(i,k) 784 D_q1(i,k) = q1(i,k-1)-q1(i,k) 785 D_q2(i,k) = q2(i,k-1)-q2(i,k) 786 D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k) 787 ENDIF 788 ENDDO 789 ENDDO 790 791 DO i=1,klon 792 omgbdth(i,1) = 0. 793 omgbdq(i,1) = 0. 794 ENDDO 795 796 DO k= 2,klev 797 DO i = 1,klon 798 IF(k .LE. kupper(i)+1) THEN ! loop on interfaces 799 omgbdth(i,k) = omgb(i,k)*( dth(i,k-1) - dth(i,k)) 800 omgbdq(i,k) = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k)) 801 ENDIF 802 ENDDO 803 ENDDO 804 c 805 c----------------------------------------------------------------- 806 DO k= 1,klev 807 DO i = 1,klon 808 IF(k .LE. kupper(i)-1) THEN 809 c----------------------------------------------------------------- 810 c 811 c Compute redistribution (advective) term 812 c 813 d_deltatw(i,k) = 814 $ dtimesub/(Ph(i,k)-Ph(i,k+1))*( 815 $ RRd1*omg(i,k )*sigmaw(i) *D_Th1(i,k) 816 $ -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) 817 $ -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)* 818 $ omgbdth(i,k+1))*ppi(i,k) 819 c print*,'d_deltatw=',d_deltatw(i,k) 820 c 821 d_deltaqw(i,k) = 822 $ dtimesub/(Ph(i,k)-Ph(i,k+1))*( 823 $ RRd1*omg(i,k )*sigmaw(i) *D_q1(i,k) 824 $ -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) 825 $ -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)* 826 $ omgbdq(i,k+1)) 827 c print*,'d_deltaqw=',d_deltaqw(i,k) 828 c 829 c and increment large scale tendencies 830 c 831 dtls(i,k) = dtls(i,k) + 832 $ dtimesub*( 833 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_Th1(i,k) 834 $ -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) ) 835 $ /(Ph(i,k)-Ph(i,k+1)) 836 $ -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k) 837 $ )*ppi(i,k) 838 c print*,'dtls=',dtls(i,k) 839 c 840 dqls(i,k) = dqls(i,k) + 841 $ dtimesub*( 842 $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_q1(i,k) 843 $ -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) ) 844 $ /(Ph(i,k)-Ph(i,k+1)) 845 $ -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k) 846 $ ) 847 c print*,'dqls=',dqls(k) 848 ENDIF 849 c------------------------------------------------------------------- 850 ENDDO 851 ENDDO 852 c------------------------------------------------------------------ 853 C 854 C Increment state variables 855 856 DO k= 1,klev 857 DO i = 1,klon 858 IF(k .LE. kupper(i)-1) THEN 859 c 860 c Coefficient de répartition 861 862 Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i)) 863 $ -ph(i,1)) 864 Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)- 865 $ ph(i,kupper(i))) 866 867 868 c Reintroduce compensating subsidence term. 869 870 c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw 871 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) 872 c . /(1-sigmaw) 873 c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw 874 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) 875 c . /(1-sigmaw) 876 c 877 c dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw 878 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) 879 c . /(1-sigmaw) 880 c dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw 881 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) 882 c . /(1-sigmaw) 883 884 dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i))) 885 dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i))) 886 c print*,'dtKE=',dtKE(k) 887 c print*,'dqKE=',dqKE(k) 888 c 889 dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i))) 890 dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i))) 891 c 892 spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ 893 $ sigmaw(i) 894 895 896 c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei 897 898 d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)* 899 $ dtimesub 900 ff(i)=d_deltatw(i,k)/dtimesub 901 902 c Sans GW 903 c 904 c deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) 905 c 906 c GW formule 1 907 c 908 c deltatw(k) = deltatw(k)+dtimesub* 909 c $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) 910 c 911 c GW formule 2 912 913 IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN 914 deltatw(i,k) = deltatw(i,k)+dtimesub* 915 $ (ff(i)+dtKE(i,k)+dtPBL(i,k) 916 $ - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k)) 917 ELSE 918 deltatw(i,k) = deltatw(i,k)+1/Tgw(i,k)*(1-exp(-dtimesub* 919 $ Tgw(i,k)))* 920 $ (ff(i)+dtKE(i,k)+dtPBL(i,k) 921 $ - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k)) 922 ENDIF 923 924 dth(i,k) = deltatw(i,k)/ppi(i,k) 925 926 gg(i)=d_deltaqw(i,k)/dtimesub 927 928 deltaqw(i,k) = deltaqw(i,k) + 929 $ dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k) - spread(i,k)* 930 $ deltaqw(i,k)) 931 932 d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) 933 d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) 934 ENDIF 935 ENDDO 936 ENDDO 937 938 C And update large scale variables 939 cIM 060208 manque DO i + remplace DO k=1,kupper(i) 940 cIM 060208 DO k = 1,kupper(i) 941 DO k= 1,klev 942 DO i = 1,klon 943 IF(k .LE. kupper(i)) THEN 944 te(i,k) = te0(i,k) + dtls(i,k) 945 qe(i,k) = qe0(i,k) + dqls(i,k) 946 the(i,k) = te(i,k)/ppi(i,k) 947 ENDIF 948 ENDDO 949 ENDDO 950 c 951 C 952 c Determine Ptop from buoyancy integral 953 c --------------------------------------- 954 c 955 c- 1/ Pressure of the level where dth changes sign. 956 c 957 DO i=1,klon 958 Ptop_provis(i)=ph(i,1) 959 ENDDO 960 c 961 DO k= 2,klev 962 DO i=1,klon 963 IF (Ptop_provis(i) .EQ. ph(i,1) .AND. 964 $ dth(i,k) .GT. -delta_t_min .and. 965 $ dth(i,k-1).LT. -delta_t_min) THEN 966 Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) 967 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k) 968 $ - dth(i,k-1)) 969 ENDIF 970 ENDDO 971 ENDDO 972 c 973 c- 2/ dth integral 974 c 975 DO i=1,klon 976 sum_dth(i) = 0. 977 dthmin(i) = -delta_t_min 978 z(i) = 0. 979 ENDDO 980 981 DO k = 1,klev 982 DO i=1,klon 983 dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg) 984 IF (dz(i) .gt. 0) THEN 985 z(i) = z(i)+dz(i) 986 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 987 dthmin(i) = amin1(dthmin(i),dth(i,k)) 988 ENDIF 989 ENDDO 990 ENDDO 991 c 992 c- 3/ height of triangle with area= sum_dth and base = dthmin 993 994 DO i=1,klon 995 hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5) 996 hw(i) = amax1(hwmin,hw(i)) 997 ENDDO 998 c 999 c- 4/ now, get Ptop 1000 c 1001 DO i=1,klon 1002 ktop(i) = 0 1003 z(i)=0. 1004 ENDDO 1005 c 1006 DO k = 1,klev 1007 DO i=1,klon 1008 dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i)) 1009 IF (dz(i) .gt. 0) THEN 1010 z(i) = z(i)+dz(i) 1011 Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i) 1012 ktop(i) = k 1013 ENDIF 1014 ENDDO 1015 ENDDO 1016 c 1017 c 4.5/Correct ktop and ptop 1018 c 1019 DO i=1,klon 1020 Ptop_new(i)=ptop(i) 1021 ENDDO 1022 c 1023 DO k= klev,2,-1 1024 DO i=1,klon 1025 cIM v3JYG; IF (k .GE. ktop(i) 1026 IF (k .LE. ktop(i) .AND. 1027 $ ptop_new(i) .EQ. ptop(i) .AND. 1028 $ dth(i,k) .GT. -delta_t_min .and. 1029 $ dth(i,k-1).LT. -delta_t_min) THEN 1030 Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) 1031 $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k) 1032 $ - dth(i,k-1)) 1033 ENDIF 1034 ENDDO 1035 ENDDO 1036 c 1037 c 1038 DO i=1,klon 1039 ptop(i) = ptop_new(i) 1040 ENDDO 1041 1042 DO k=klev,1,-1 1043 DO i=1,klon 1044 IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k 1045 ENDDO 1046 ENDDO 1047 c 1048 c 5/ Set deltatw & deltaqw to 0 above kupper 1049 c 1050 DO k = 1,klev 1051 DO i=1,klon 1052 IF (k .GE. kupper(i)) THEN 1053 deltatw(i,k) = 0. 1054 deltaqw(i,k) = 0. 1055 ENDIF 1056 ENDDO 1057 ENDDO 1058 c 1059 C 1060 ENDDO ! end sub-timestep loop 1061 C 1062 C ----------------------------------------------------------------- 1063 c Get back to tendencies per second 1064 c 1065 DO k = 1,klev 1066 DO i=1,klon 1067 IF (k .LE. kupper(i)-1) THEN 1068 dtls(i,k) = dtls(i,k)/dtime 1069 dqls(i,k) = dqls(i,k)/dtime 1070 d_deltatw2(i,k)=d_deltatw2(i,k)/dtime 1071 d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime 1072 d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime 1073 ENDIF 1074 ENDDO 1075 ENDDO 1076 c 1077 c 1078 c---------------------------------------------------------- 1079 c Determine wake final state; recompute wape, cstar, ktop; 1080 c filter out bad wakes. 1081 c---------------------------------------------------------- 1082 c 1083 C 2.1 - Undisturbed area and Wake integrals 1084 C --------------------------------------------------------- 1085 1086 DO i=1,klon 1087 z(i) = 0. 1088 sum_thu(i) = 0. 1089 sum_tu(i) = 0. 1090 sum_qu(i) = 0. 1091 sum_thvu(i) = 0. 1092 sum_dth(i) = 0. 1093 sum_dq(i) = 0. 1094 sum_rho(i) = 0. 1095 sum_dtdwn(i) = 0. 1096 sum_dqdwn(i) = 0. 1097 1098 av_thu(i) = 0. 1099 av_tu(i) =0. 1100 av_qu(i) =0. 1101 av_thvu(i) = 0. 1102 av_dth(i) = 0. 1103 av_dq(i) = 0. 1104 av_rho(i) =0. 1105 av_dtdwn(i) =0. 1106 av_dqdwn(i) = 0. 1107 ENDDO 1108 C Potential temperatures and humidity 1109 c---------------------------------------------------------- 1110 1111 DO k =1,klev 1112 DO i=1,klon 1113 rho(i,k) = p(i,k)/(rd*te(i,k)) 1114 IF(k .eq. 1) THEN 1115 rhoh(i,k) = ph(i,k)/(rd*te(i,k)) 1116 zhh(i,k)=0 1117 ELSE 1118 rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1))) 1119 zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1) 1120 ENDIF 1121 the(i,k) = te(i,k)/ppi(i,k) 1122 thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k) 1123 tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i) 1124 qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i) 1125 rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k))) 1126 dth(i,k) = deltatw(i,k)/ppi(i,k) 1127 ENDDO 1128 ENDDO 1129 1130 C Integrals (and wake top level number) 1131 C ----------------------------------------------------------- 1132 1133 C Initialize sum_thvu to 1st level virt. pot. temp. 1134 1135 DO i=1,klon 1136 z(i) = 1. 1137 dz(i) = 1. 1138 sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) 1139 sum_dth(i) = 0. 1140 ENDDO 1141 1142 DO k = 1,klev 1143 DO i=1,klon 1144 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) 1145 IF (dz(i) .GT. 0) THEN 1146 z(i) = z(i)+dz(i) 1147 sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i) 1148 sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i) 1149 sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i) 1150 sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i) 1151 sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) 1152 sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i) 1153 sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i) 1154 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i) 1155 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i) 1156 ENDIF 1157 ENDDO 1158 ENDDO 1159 c 1160 DO i=1,klon 1161 hw0(i) = z(i) 1162 ENDDO 1163 c 1164 C 2.1 - WAPE and mean forcing computation 1165 C------------------------------------------------------------- 1166 1167 C Means 1168 1169 DO i=1, klon 1170 av_thu(i) = sum_thu(i)/hw0(i) 1171 av_tu(i) = sum_tu(i)/hw0(i) 1172 av_qu(i) = sum_qu(i)/hw0(i) 1173 av_thvu(i) = sum_thvu(i)/hw0(i) 1174 av_dth(i) = sum_dth(i)/hw0(i) 1175 av_dq(i) = sum_dq(i)/hw0(i) 1176 av_rho(i) = sum_rho(i)/hw0(i) 1177 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 1178 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 1179 1180 wape2(i) = - rg*hw0(i)*(av_dth(i) 1181 $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+ 1182 $ av_dth(i)*av_dq(i) ))/av_thvu(i) 1183 ENDDO 1184 1185 C 2.2 Prognostic variable update 1186 C ------------------------------------------------------------ 1187 1188 C Filter out bad wakes 1189 c 1190 DO k = 1,klev 1191 DO i=1,klon 1192 IF ( wape2(i) .LT. 0.) THEN 1193 deltatw(i,k) = 0. 1194 deltaqw(i,k) = 0. 1195 dth(i,k) = 0. 1196 ENDIF 1197 ENDDO 1198 ENDDO 1199 c 1200 1201 DO i=1, klon 1202 IF ( wape2(i) .LT. 0.) THEN 1203 wape2(i) = 0. 1204 Cstar2(i) = 0. 1205 hw(i) = hwmin 1206 sigmaw(i) = amax1(sigmad,sigd_con(i)) 1207 fip(i) = 0. 1208 gwake(i) = .FALSE. 1209 ELSE 1210 if(prt_level.ge.10) print*,'wape2>0' 1211 Cstar2(i) = stark*sqrt(2.*wape2(i)) 1212 gwake(i) = .TRUE. 1213 ENDIF 1214 ENDDO 1215 c 1216 DO i=1, klon 1217 ktopw(i) = ktop(i) 1218 ENDDO 1219 c 1220 DO i=1, klon 1221 IF (ktopw(i) .gt. 0 .and. gwake(i)) then 1222 1223 Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) 1224 ccc heff = 600. 1225 C Utilisation de la hauteur hw 1226 cc heff = 0.7*hw 1227 heff(i) = hw(i) 1228 1229 FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2* 1230 $ sqrt(sigmaw(i)*wdens*3.14) 1231 FIP(i) = alpk * FIP(i) 1232 Cjyg2 1233 ELSE 1234 FIP(i) = 0. 1235 ENDIF 1236 ENDDO 1237 c 1238 C Limitation de sigmaw 1239 c 1240 C sécurité : si le wake occuppe plus de 90 % de la surface de la maille, 1241 C alors il disparait en se mélangeant à la partie undisturbed 1242 c 1243 DO k = 1,klev 1244 DO i=1, klon 1245 IF ((sigmaw(i).GT.0.9).or. 1246 $ ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN 1247 ccc IF (sigmaw(i).GT.0.9) THEN 1248 dtls(i,k) = 0. 1249 dqls(i,k) = 0. 1250 deltatw(i,k) = 0. 1251 deltaqw(i,k) = 0. 1252 ENDIF 1253 ENDDO 1254 ENDDO 1255 c 1256 DO i=1, klon 1257 IF ((sigmaw(i).GT.0.9).or. 1258 $ ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN 1259 ccc IF (sigmaw(i).GT.0.9) THEN 1260 wape(i) = 0. 1261 hw(i) = hwmin 1262 sigmaw(i) = sigmad 1263 fip(i) = 0. 1264 ELSE 1265 wape(i) = wape2(i) 1266 ENDIF 1267 ENDDO 1268 c 1269 c 1270 RETURN 1271 END 1272 Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con 1273 : ,te0,qe0,omgb 1274 : ,dtdwn,dqdwn,amdwn,amup,dta,dqa 1275 : ,wdtPBL,wdqPBL,udtPBL,udqPBL 1276 o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl 1277 o ,dtls,dqls 1278 o ,ktopw,omgbdth,dp_omgb,wdens 1279 o ,tu,qu 1280 o ,dtKE,dqKE 1281 o ,dtPBL,dqPBL 1282 o ,omg,dp_deltomg,spread 1283 o ,Cstar,d_deltat_gw 1284 o ,d_deltatw2,d_deltaqw2) 1285 1286 *************************************************************** 1287 * * 1288 * WAKE * 1289 * retour a un Pupper fixe * 1290 * * 1291 * written by : GRANDPEIX Jean-Yves 09/03/2000 * 1292 * modified by : ROEHRIG Romain 01/29/2007 * 1293 *************************************************************** 1294 c 1295 USE dimphy 1296 IMPLICIT none 1297 c============================================================================ 1298 C 1299 C 1300 C But : Decrire le comportement des poches froides apparaissant dans les 1301 C grands systemes convectifs, et fournir l'energie disponible pour 1302 C le declenchement de nouvelles colonnes convectives. 1303 C 1304 C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area 1305 C deltaqw : ecart d'humidite wake-undisturbed area 1306 C sigmaw : fraction d'aire occupee par la poche. 1307 C 1308 C Variable de sortie : 1309 c 1310 c wape : WAke Potential Energy 1311 c fip : Front Incident Power (W/m2) - ALP 1312 c gfl : Gust Front Length per unit area (m-1) 1313 C dtls : large scale temperature tendency due to wake 1314 C dqls : large scale humidity tendency due to wake 1315 C hw : hauteur de la poche 1316 C dp_omgb : vertical gradient of large scale omega 1317 C omgbdth: flux of Delta_Theta transported by LS omega 1318 C dtKE : differential heating (wake - unpertubed) 1319 C dqKE : differential moistening (wake - unpertubed) 1320 C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) 1321 C dp_deltomg : vertical gradient of omg (s-1) 1322 C spread : spreading term in dt_wake and dq_wake 1323 C deltatw : updated temperature difference (T_w-T_u). 1324 C deltaqw : updated humidity difference (q_w-q_u). 1325 C sigmaw : updated wake fractional area. 1326 C d_deltat_gw : delta T tendency due to GW 1327 c 1328 C Variables d'entree : 1329 c 1330 c aire : aire de la maille 1331 c te0 : temperature dans l'environnement (K) 1332 C qe0 : humidite dans l'environnement (kg/kg) 1333 C omgb : vitesse verticale moyenne sur la maille (Pa/s) 1334 C dtdwn: source de chaleur due aux descentes (K/s) 1335 C dqdwn: source d'humidite due aux descentes (kg/kg/s) 1336 C dta : source de chaleur due courants satures et detrain (K/s) 1337 C dqa : source d'humidite due aux courants satures et detra (kg/kg/s) 1338 C amdwn: flux de masse total des descentes, par unite de 1339 C surface de la maille (kg/m2/s) 1340 C amup : flux de masse total des ascendances, par unite de 1341 C surface de la maille (kg/m2/s) 1342 C p : pressions aux milieux des couches (Pa) 1343 C ph : pressions aux interfaces (Pa) 1344 C ppi : (p/p_0)**kapa (adim) 1345 C dtime: increment temporel (s) 1346 c 1347 C Variables internes : 1348 c 1349 c rhow : masse volumique de la poche froide 1350 C rho : environment density at P levels 1351 C rhoh : environment density at Ph levels 1352 C te : environment temperature | may change within 1353 C qe : environment humidity | sub-time-stepping 1354 C the : environment potential temperature 1355 C thu : potential temperature in undisturbed area 1356 C tu : temperature in undisturbed area 1357 C qu : humidity in undisturbed area 1358 C dp_omgb: vertical gradient og LS omega 1359 C omgbw : wake average vertical omega 1360 C dp_omgbw: vertical gradient of omgbw 1361 C omgbdq : flux of Delta_q transported by LS omega 1362 C dth : potential temperature diff. wake-undist. 1363 C th1 : first pot. temp. for vertical advection (=thu) 1364 C th2 : second pot. temp. for vertical advection (=thw) 1365 C q1 : first humidity for vertical advection 1366 C q2 : second humidity for vertical advection 1367 C d_deltatw : terme de redistribution pour deltatw 1368 C d_deltaqw : terme de redistribution pour deltaqw 1369 C deltatw0 : deltatw initial 1370 C deltaqw0 : deltaqw initial 1371 C hw0 : hw initial 1372 C sigmaw0: sigmaw initial 1373 C amflux : horizontal mass flux through wake boundary 1374 C wdens : number of wakes per unit area (3D) or per 1375 C unit length (2D) 1376 C Tgw : 1 sur la période de onde de gravité 1377 c Cgw : vitesse de propagation de onde de gravité 1378 c LL : distance entre 2 poches 1379 1380 c------------------------------------------------------------------------- 1381 c Déclaration de variables 1382 c------------------------------------------------------------------------- 1383 1384 #include "dimensions.h" 114 1385 cccc#include "dimphy.h" 115 1386 #include "YOMCST.h"
Note: See TracChangeset
for help on using the changeset viewer.
