| 1 | c $Header$ |
|---|
| 2 | SUBROUTINE CONVECT |
|---|
| 3 | * (DTIME, T1, R1, RS, U, V, TRA, P, PH, |
|---|
| 4 | * ND, NDP1, NL, NTRA, DELT, IFLAG, |
|---|
| 5 | * FT, FR, FU, FV, FTRA, PRECIP, |
|---|
| 6 | * icb,inb, upwd,dnwd,dnwd0,SIG, W0,Mike,Mke, |
|---|
| 7 | * Ma,MENTS,QENTS,TPS,TLS,SIGIJ,CAPE,TVP,PBASE,BUOYBASE, |
|---|
| 8 | * DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR) |
|---|
| 9 | C |
|---|
| 10 | C *** THE PARAMETER NA SHOULD IN GENERAL EQUAL ND *** |
|---|
| 11 | C |
|---|
| 12 | c################################################################# |
|---|
| 13 | cFleur Introduction des traceurs dans convect3 le 6 juin 200 |
|---|
| 14 | c################################################################# |
|---|
| 15 | #include "dimensions.h" |
|---|
| 16 | #include "dimphy.h" |
|---|
| 17 | PARAMETER (NA=60) |
|---|
| 18 | |
|---|
| 19 | integer NTRAC |
|---|
| 20 | PARAMETER (NTRAC=nqmx-2) |
|---|
| 21 | |
|---|
| 22 | INTEGER NENT(NA) |
|---|
| 23 | REAL T1(ND),R1(ND),RS(ND),U(ND),V(ND),TRA(ND,NTRA) |
|---|
| 24 | REAL P(ND),PH(NDP1) |
|---|
| 25 | REAL FT(ND),FR(ND),FU(ND),FV(ND),FTRA(ND,NTRA) |
|---|
| 26 | REAL SIG(ND),W0(ND) |
|---|
| 27 | REAL UENT(NA,NA),VENT(NA,NA),TRAENT(NA,NA,NTRAC),TRATM(NA) |
|---|
| 28 | REAL UP(NA),VP(NA),TRAP(NA,NTRAC) |
|---|
| 29 | REAL M(NA),MP(NA),MENT(NA,NA),QENT(NA,NA),ELIJ(NA,NA) |
|---|
| 30 | REAL SIJ(NA,NA),TVP(NA),TV(NA),WATER(NA) |
|---|
| 31 | REAL RP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA) |
|---|
| 32 | REAL SIGP(NA),B(NA),TP(NA),CPN(NA) |
|---|
| 33 | REAL LV(NA),LVCP(NA),H(NA),HP(NA),GZ(NA) |
|---|
| 34 | REAL T(NA),RR(NA) |
|---|
| 35 | C |
|---|
| 36 | REAL BUOY(NA) ! Lifted parcel buoyancy |
|---|
| 37 | REAL DTVPDT1(ND),DTVPDQ1(ND) ! Derivatives of parcel virtual |
|---|
| 38 | C temperature wrt T1 and Q1 |
|---|
| 39 | REAL CLW_NEW(NA),QI(NA) |
|---|
| 40 | C |
|---|
| 41 | LOGICAL ICE_CONV |
|---|
| 42 | PARAMETER (ICE_CONV=.TRUE.) |
|---|
| 43 | |
|---|
| 44 | cccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 45 | c declaration des variables a sortir |
|---|
| 46 | ccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 47 | real Mke(nd) |
|---|
| 48 | real Mike(nd) |
|---|
| 49 | real Ma(nd) |
|---|
| 50 | real TPS(ND) !temperature dans les ascendances non diluees |
|---|
| 51 | real TLS(ND) !temperature potentielle |
|---|
| 52 | real MENTS(nd,nd) |
|---|
| 53 | real QENTS(nd,nd) |
|---|
| 54 | REAL SIGIJ(KLEV,KLEV) |
|---|
| 55 | REAL PBASE ! pressure at the cloud base level |
|---|
| 56 | REAL BUOYBASE ! buoyancy at the cloud base level |
|---|
| 57 | cccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 58 | |
|---|
| 59 | |
|---|
| 60 | |
|---|
| 61 | c |
|---|
| 62 | real dnwd0(nd) ! precipitation driven unsaturated downdraft flux |
|---|
| 63 | real dnwd(nd), dn1 ! in-cloud saturated downdraft mass flux |
|---|
| 64 | real upwd(nd), up1 ! in-cloud saturated updraft mass flux |
|---|
| 65 | C |
|---|
| 66 | C *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *** |
|---|
| 67 | C *** THESE SHOULD BE CONSISTENT WITH *** |
|---|
| 68 | C *** THOSE USED IN CALLING PROGRAM *** |
|---|
| 69 | C *** NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT *** |
|---|
| 70 | C |
|---|
| 71 | c sb CPD=1005.7 |
|---|
| 72 | c sb CPV=1870.0 |
|---|
| 73 | c sb CL=4190.0 |
|---|
| 74 | c sb CPVMCL=CL-CPV |
|---|
| 75 | c sb RV=461.5 |
|---|
| 76 | c sb RD=287.04 |
|---|
| 77 | c sb EPS=RD/RV |
|---|
| 78 | c sb ALV0=2.501E6 |
|---|
| 79 | ccccccccccccccccccccccc |
|---|
| 80 | c constantes coherentes avec le modele du Centre Europeen |
|---|
| 81 | c sb RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644 |
|---|
| 82 | c sb RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153 |
|---|
| 83 | c sb CPD = 3.5 * RD |
|---|
| 84 | c sb CPV = 4.0 * RV |
|---|
| 85 | c sb CL = 4218.0 |
|---|
| 86 | c sb CPVMCL=CL-CPV |
|---|
| 87 | c sb EPS=RD/RV |
|---|
| 88 | c sb ALV0=2.5008E+06 |
|---|
| 89 | cccccccccccccccccccccc |
|---|
| 90 | c on utilise les constantes thermo du Centre Europeen: (SB) |
|---|
| 91 | c |
|---|
| 92 | #include "YOMCST.h" |
|---|
| 93 | c |
|---|
| 94 | CPD = RCPD |
|---|
| 95 | CPV = RCPV |
|---|
| 96 | CL = RCW |
|---|
| 97 | CPVMCL = CL-CPV |
|---|
| 98 | EPS = RD/RV |
|---|
| 99 | ALV0 = RLVTT |
|---|
| 100 | c |
|---|
| 101 | NK = 1 ! origin level of the lifted parcel |
|---|
| 102 | c |
|---|
| 103 | cccccccccccccccccccccc |
|---|
| 104 | C |
|---|
| 105 | C *** INITIALIZE OUTPUT ARRAYS AND PARAMETERS *** |
|---|
| 106 | C |
|---|
| 107 | DO 5 I=1,ND |
|---|
| 108 | FT(I)=0.0 |
|---|
| 109 | FR(I)=0.0 |
|---|
| 110 | FU(I)=0.0 |
|---|
| 111 | FV(I)=0.0 |
|---|
| 112 | DO 4 J=1,NTRA |
|---|
| 113 | FTRA(I,J)=0.0 |
|---|
| 114 | 4 CONTINUE |
|---|
| 115 | T(I)=T1(I) |
|---|
| 116 | RR(I)=R1(I) |
|---|
| 117 | 5 CONTINUE |
|---|
| 118 | DO 7 I=1,NL |
|---|
| 119 | RDCP=(RD*(1.-RR(I))+RR(I)*RV)/ (CPD*(1.-RR(I))+RR(I)*CPV) |
|---|
| 120 | TH(I)=T(I)*(1000.0/P(I))**RDCP |
|---|
| 121 | 7 CONTINUE |
|---|
| 122 | C |
|---|
| 123 | ************************************************************* |
|---|
| 124 | ** CALCUL DES TEMPERATURES POTENTIELLES A SORTIR |
|---|
| 125 | ************************************************************* |
|---|
| 126 | do i=1,ND |
|---|
| 127 | RDCP=(RD*(1.-RR(I))+RR(I)*RV)/ (CPD*(1.-RR(I))+RR(I)*CPV) |
|---|
| 128 | |
|---|
| 129 | TLS(i)=T(I)*(1000.0/P(I))**RDCP |
|---|
| 130 | enddo |
|---|
| 131 | |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | |
|---|
| 135 | ************************************************************ |
|---|
| 136 | |
|---|
| 137 | |
|---|
| 138 | PRECIP=0.0 |
|---|
| 139 | IFLAG=1 |
|---|
| 140 | C |
|---|
| 141 | C *** SPECIFY PARAMETERS *** |
|---|
| 142 | C *** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *** |
|---|
| 143 | C *** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *** |
|---|
| 144 | C *** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *** |
|---|
| 145 | C *** EFFICIENCY IS ASSUMED TO BE UNITY *** |
|---|
| 146 | C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** |
|---|
| 147 | C *** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** |
|---|
| 148 | C *** OF CLOUD *** |
|---|
| 149 | C *** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *** |
|---|
| 150 | C *** APPROACH TO QUASI-EQUILIBRIUM *** |
|---|
| 151 | C *** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *** |
|---|
| 152 | C *** (BETA MUST BE LESS THAN OR EQUAL TO 1) *** |
|---|
| 153 | C *** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *** |
|---|
| 154 | C *** APPROACH TO QUASI-EQUILIBRIUM *** |
|---|
| 155 | C *** IT MUST BE LESS THAN 0 *** |
|---|
| 156 | C |
|---|
| 157 | PBCRIT=150.0 |
|---|
| 158 | PTCRIT=500.0 |
|---|
| 159 | SIGD=0.01 |
|---|
| 160 | SPFAC=0.15 |
|---|
| 161 | c sb: |
|---|
| 162 | EPMAX=0.993 ! precip efficiency less than unity |
|---|
| 163 | c EPMAX=1. ! precip efficiency less than unity |
|---|
| 164 | C |
|---|
| 165 | Cjyg |
|---|
| 166 | CCC BETA=0.96 |
|---|
| 167 | C Beta is now expressed as a function of the characteristic time |
|---|
| 168 | C of the convective process. |
|---|
| 169 | CCC Old value : TAU = 15000. !(for dtime = 600.s) |
|---|
| 170 | CCC Other value (inducing little change) :TAU = 8000. |
|---|
| 171 | TAU = 8000. |
|---|
| 172 | BETA = 1.-DTIME/TAU |
|---|
| 173 | Cjyg |
|---|
| 174 | CCC ALPHA=1.0 |
|---|
| 175 | ALPHA=1.5E-3*DTIME/TAU |
|---|
| 176 | C Increase alpha in order to compensate W decrease |
|---|
| 177 | ALPHA = ALPHA*1.5 |
|---|
| 178 | C |
|---|
| 179 | Cjyg (voir CONVECT 3) |
|---|
| 180 | CCC DTCRIT=-0.2 |
|---|
| 181 | DTCRIT=-2. |
|---|
| 182 | Cgf&jyg |
|---|
| 183 | CCC DT pour l'overshoot. |
|---|
| 184 | DTOVSH = -0.2 |
|---|
| 185 | |
|---|
| 186 | C |
|---|
| 187 | C *** INCREMENT THE COUNTER *** |
|---|
| 188 | C |
|---|
| 189 | SIG(ND)=SIG(ND)+1.0 |
|---|
| 190 | SIG(ND)=AMIN1(SIG(ND),12.1) |
|---|
| 191 | C |
|---|
| 192 | C *** IF NOPT IS AN INTEGER OTHER THAN 0, CONVECT *** |
|---|
| 193 | C *** RETURNS ARRAYS T AND R THAT MAY HAVE BEEN *** |
|---|
| 194 | C *** ALTERED BY DRY ADIABATIC ADJUSTMENT; OTHERWISE *** |
|---|
| 195 | C *** THE RETURNED ARRAYS ARE UNALTERED. *** |
|---|
| 196 | C |
|---|
| 197 | NOPT=0 |
|---|
| 198 | C |
|---|
| 199 | C *** PERFORM DRY ADIABATIC ADJUSTMENT *** |
|---|
| 200 | C |
|---|
| 201 | C *** DO NOT BYPASS THIS EVEN IF THE CALLING PROGRAM HAS A *** |
|---|
| 202 | C *** BOUNDARY LAYER SCHEME !!! *** |
|---|
| 203 | C |
|---|
| 204 | c test sb DO 30 I=NL-1,1,-1 |
|---|
| 205 | c test sb JN=0 |
|---|
| 206 | c test sb DO 10 J=I+1,NL |
|---|
| 207 | c test sb 10 IF(TH(J).LT.TH(I))JN=J |
|---|
| 208 | c test sb IF(JN.EQ.0)GOTO 30 |
|---|
| 209 | c test sb AHM=0.0 |
|---|
| 210 | c test sb RM=0.0 |
|---|
| 211 | c test sb UM=0.0 |
|---|
| 212 | c test sb VM=0.0 |
|---|
| 213 | c test sb DO K=1,NTRA |
|---|
| 214 | c test sb TRATM(K)=0.0 |
|---|
| 215 | c test sb END DO |
|---|
| 216 | c test sb DO 15 J=I,JN |
|---|
| 217 | c test sb AHM=AHM+(CPD*(1.-RR(J))+RR(J)*CPV)*T(J)*(PH(J)-PH(J+1)) |
|---|
| 218 | c test sb RM=RM+RR(J)*(PH(J)-PH(J+1)) |
|---|
| 219 | c test sb UM=UM+U(J)*(PH(J)-PH(J+1)) |
|---|
| 220 | c test sb VM=VM+V(J)*(PH(J)-PH(J+1)) |
|---|
| 221 | c test sb DO K=1,NTRA |
|---|
| 222 | c test sb TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1)) |
|---|
| 223 | c test sb END DO |
|---|
| 224 | c test sb 15 CONTINUE |
|---|
| 225 | c test sb DPHINV=1./(PH(I)-PH(JN+1)) |
|---|
| 226 | c test sb RM=RM*DPHINV |
|---|
| 227 | c test sb UM=UM*DPHINV |
|---|
| 228 | c test sb VM=VM*DPHINV |
|---|
| 229 | c test sb DO K=1,NTRA |
|---|
| 230 | c test sb TRATM(K)=TRATM(K)*DPHINV |
|---|
| 231 | c test sb END DO |
|---|
| 232 | c test sb A2=0.0 |
|---|
| 233 | c test sb DO 20 J=I,JN |
|---|
| 234 | c test sb RR(J)=RM |
|---|
| 235 | c test sb U(J)=UM |
|---|
| 236 | c test sb V(J)=VM |
|---|
| 237 | c test sb DO K=1,NTRA |
|---|
| 238 | c test sb TRA(J,K)=TRATM(K) |
|---|
| 239 | c test sb END DO |
|---|
| 240 | c test sb RDCP=(RD*(1.-RR(J))+RR(J)*RV)/ (CPD*(1.-RR(J))+RR(J)*CPV) |
|---|
| 241 | c test sb X=(0.001*P(J))**RDCP |
|---|
| 242 | c test sb T(J)=X |
|---|
| 243 | c test sb A2=A2+(CPD*(1.-RR(J))+RR(J)*CPV)*X*(PH(J)-PH(J+1)) |
|---|
| 244 | c test sb 20 CONTINUE |
|---|
| 245 | c test sb DO 25 J=I,JN |
|---|
| 246 | c test sb TH(J)=AHM/A2 |
|---|
| 247 | c test sb T(J)=T(J)*TH(J) |
|---|
| 248 | c test sb 25 CONTINUE |
|---|
| 249 | c test sb 30 CONTINUE |
|---|
| 250 | C |
|---|
| 251 | C *** RESET INPUT ARRAYS IF NOPT .NE. 0 *** |
|---|
| 252 | C |
|---|
| 253 | IF (NOPT.NE.0)THEN |
|---|
| 254 | DO 35 I=1,ND |
|---|
| 255 | T1(I)=T(I) |
|---|
| 256 | R1(I)=RR(I) |
|---|
| 257 | 35 CONTINUE |
|---|
| 258 | END IF |
|---|
| 259 | C |
|---|
| 260 | C *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY |
|---|
| 261 | C |
|---|
| 262 | GZ(1)=0.0 |
|---|
| 263 | CPN(1)=CPD*(1.-RR(1))+RR(1)*CPV |
|---|
| 264 | H(1)=T(1)*CPN(1) |
|---|
| 265 | DO 40 I=2,NL |
|---|
| 266 | TVX=T(I)*(1.+RR(I)/EPS-RR(I)) |
|---|
| 267 | TVY=T(I-1)*(1.+RR(I-1)/EPS-RR(I-1)) |
|---|
| 268 | GZ(I)=GZ(I-1)+0.5*RD*(TVX+TVY)*(P(I-1)-P(I))/PH(I) |
|---|
| 269 | CPN(I)=CPD*(1.-RR(I))+CPV*RR(I) |
|---|
| 270 | H(I)=T(I)*CPN(I)+GZ(I) |
|---|
| 271 | 40 CONTINUE |
|---|
| 272 | C |
|---|
| 273 | C *** CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT LOWEST MODEL LEVEL *** |
|---|
| 274 | C *** (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980) *** |
|---|
| 275 | C |
|---|
| 276 | IF (T(1).LT.250.0.OR.RR(1).LE.0.0)THEN |
|---|
| 277 | IFLAG=0 |
|---|
| 278 | c sb3d print*,'je suis passe par 366' |
|---|
| 279 | RETURN |
|---|
| 280 | END IF |
|---|
| 281 | |
|---|
| 282 | cjyg1 Utilisation de la subroutine CLIFT |
|---|
| 283 | CC RH=RR(1)/RS(1) |
|---|
| 284 | CC CHI=T(1)/(1669.0-122.0*RH-T(1)) |
|---|
| 285 | CC PLCL=P(1)*(RH**CHI) |
|---|
| 286 | CALL CLIFT(P(1),T(1),RR(1),RS(1),PLCL,DPLCLDT,DPLCLDR) |
|---|
| 287 | cjyg2 |
|---|
| 288 | c sb3d PRINT *,' em_plcl,p1,t1,r1,rs1,rh ' |
|---|
| 289 | c sb3d $ ,PLCL,P(1),T(1),RR(1),RS(1),RH |
|---|
| 290 | c |
|---|
| 291 | IF (PLCL.LT.200.0.OR.PLCL.GE.2000.0)THEN |
|---|
| 292 | IFLAG=2 |
|---|
| 293 | RETURN |
|---|
| 294 | END IF |
|---|
| 295 | Cjyg1 |
|---|
| 296 | C Essais de modification de ICB |
|---|
| 297 | C |
|---|
| 298 | C *** CALCULATE FIRST LEVEL ABOVE LCL (=ICB) *** |
|---|
| 299 | C |
|---|
| 300 | CC ICB=NL-1 |
|---|
| 301 | CC DO 50 I=2,NL-1 |
|---|
| 302 | CC IF(P(I).LT.PLCL)THEN |
|---|
| 303 | CC ICB=MIN(ICB,I) ! ICB sup ou egal a 2 |
|---|
| 304 | CC END IF |
|---|
| 305 | CC 50 CONTINUE |
|---|
| 306 | CC IF(ICB.EQ.(NL-1))THEN |
|---|
| 307 | CC IFLAG=3 |
|---|
| 308 | CC RETURN |
|---|
| 309 | CC END IF |
|---|
| 310 | C |
|---|
| 311 | C *** CALCULATE LAYER CONTAINING LCL (=ICB) *** |
|---|
| 312 | C |
|---|
| 313 | ICB=NL-1 |
|---|
| 314 | c sb DO 50 I=2,NL-1 |
|---|
| 315 | DO 50 I=3,NL-1 ! modif sb pour que ICB soit sup/egal a 2 |
|---|
| 316 | C la modification consiste a comparer PLCL a PH et non a P: |
|---|
| 317 | C ICB est defini par : PH(ICB)<PLCL<PH(ICB-!) |
|---|
| 318 | IF(PH(I).LT.PLCL)THEN |
|---|
| 319 | ICB=MIN(ICB,I) |
|---|
| 320 | END IF |
|---|
| 321 | 50 CONTINUE |
|---|
| 322 | IF(ICB.EQ.(NL-1))THEN |
|---|
| 323 | IFLAG=3 |
|---|
| 324 | RETURN |
|---|
| 325 | END IF |
|---|
| 326 | ICB = ICB-1 ! ICB sup ou egal a 2 |
|---|
| 327 | Cjyg2 |
|---|
| 328 | C |
|---|
| 329 | C |
|---|
| 330 | |
|---|
| 331 | C *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL *** |
|---|
| 332 | C *** TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC *** |
|---|
| 333 | C *** LIQUID WATER CONTENT *** |
|---|
| 334 | C |
|---|
| 335 | |
|---|
| 336 | cjyg1 |
|---|
| 337 | c make sure that "Cloud base" seen by TLIFT is actually the |
|---|
| 338 | c fisrt level where adiabatic ascent is saturated |
|---|
| 339 | IF (PLCL .GT. P(ICB)) THEN |
|---|
| 340 | c sb CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,TVP,TP,CLW,ND,NL) |
|---|
| 341 | CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,NK,TVP,TP,CLW,ND,NL |
|---|
| 342 | : ,DTVPDT1,DTVPDQ1) |
|---|
| 343 | ELSE |
|---|
| 344 | c sb CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,TVP,TP,CLW,ND,NL) |
|---|
| 345 | CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,NK,TVP,TP,CLW,ND,NL |
|---|
| 346 | : ,DTVPDT1,DTVPDQ1) |
|---|
| 347 | ENDIF |
|---|
| 348 | cjyg2 |
|---|
| 349 | |
|---|
| 350 | ****************************************************************************** |
|---|
| 351 | **** SORTIE DE LA TEMPERATURE DE L ASCENDANCE NON DILUE |
|---|
| 352 | ****************************************************************************** |
|---|
| 353 | do i=1,ND |
|---|
| 354 | TPS(i)=TP(i) |
|---|
| 355 | enddo |
|---|
| 356 | |
|---|
| 357 | |
|---|
| 358 | ****************************************************************************** |
|---|
| 359 | |
|---|
| 360 | C |
|---|
| 361 | C *** SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF *** |
|---|
| 362 | C *** PRECIPITATION FALLING OUTSIDE OF CLOUD *** |
|---|
| 363 | C *** THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) *** |
|---|
| 364 | C |
|---|
| 365 | DO 55 I=1,NL |
|---|
| 366 | PDEN=PTCRIT-PBCRIT |
|---|
| 367 | c |
|---|
| 368 | cjyg |
|---|
| 369 | ccc EP(I)=(P(ICB)-P(I)-PBCRIT)/PDEN |
|---|
| 370 | c sb EP(I)=(PLCL-P(I)-PBCRIT)/PDEN |
|---|
| 371 | EP(I)=(PLCL-P(I)-PBCRIT)/PDEN * EPMAX ! sb |
|---|
| 372 | c |
|---|
| 373 | EP(I)=AMAX1(EP(I),0.0) |
|---|
| 374 | c sb EP(I)=AMIN1(EP(I),1.0) |
|---|
| 375 | EP(I)=AMIN1(EP(I),EPMAX) ! sb |
|---|
| 376 | SIGP(I)=SPFAC |
|---|
| 377 | C |
|---|
| 378 | C *** CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL *** |
|---|
| 379 | C *** VIRTUAL TEMPERATURE *** |
|---|
| 380 | C |
|---|
| 381 | TV(I)=T(I)*(1.+RR(I)/EPS-RR(I)) |
|---|
| 382 | Ccd1 |
|---|
| 383 | C . Keep all liquid water in lifted parcel (-> adiabatic CAPE) |
|---|
| 384 | C |
|---|
| 385 | ccc TVP(I)=TVP(I)-TP(I)*(RR(1)-EP(I)*CLW(I)) |
|---|
| 386 | c!!! sb TVP(I)=TVP(I)-TP(I)*RR(1) ! calcule dans tlift |
|---|
| 387 | Ccd2 |
|---|
| 388 | C |
|---|
| 389 | C *** Calculate first estimate of buoyancy |
|---|
| 390 | C |
|---|
| 391 | BUOY(I) = TVP(I) - TV(I) |
|---|
| 392 | 55 CONTINUE |
|---|
| 393 | C |
|---|
| 394 | C *** Set Cloud Base Buoyancy at (Plcl+DPbase) level buoyancy |
|---|
| 395 | C |
|---|
| 396 | DPBASE = -40. !That is 400m above LCL |
|---|
| 397 | PBASE = PLCL + DPBASE |
|---|
| 398 | TVPBASE = TVP(ICB )*(PBASE -P(ICB+1))/(P(ICB)-P(ICB+1)) |
|---|
| 399 | $ +TVP(ICB+1)*(P(ICB)-PBASE )/(P(ICB)-P(ICB+1)) |
|---|
| 400 | TVBASE = TV(ICB )*(PBASE -P(ICB+1))/(P(ICB)-P(ICB+1)) |
|---|
| 401 | $ +TV(ICB+1)*(P(ICB)-PBASE )/(P(ICB)-P(ICB+1)) |
|---|
| 402 | C |
|---|
| 403 | BUOYBASE = TVPBASE - TVBASE |
|---|
| 404 | C |
|---|
| 405 | CC Set buoyancy = BUOYBASE for all levels below BASE. |
|---|
| 406 | CC For safety, set : BUOY(ICB) = BUOYBASE |
|---|
| 407 | DO I = ICB,NL |
|---|
| 408 | IF (P(I) .GE. PBASE) THEN |
|---|
| 409 | BUOY(I) = BUOYBASE |
|---|
| 410 | ENDIF |
|---|
| 411 | ENDDO |
|---|
| 412 | BUOY(ICB) = BUOYBASE |
|---|
| 413 | C |
|---|
| 414 | c sb3d print *,'buoybase,tvp_tv,tvpbase,tvbase,pbase,plcl' |
|---|
| 415 | c sb3d $, buoybase,tvp(icb)-tv(icb),tvpbase,tvbase,pbase,plcl |
|---|
| 416 | c sb3d print *,'TVP ',(tvp(i),i=1,nl) |
|---|
| 417 | c sb3d print *,'TV ',(tv(i),i=1,nl) |
|---|
| 418 | c sb3d print *, 'P ',(p(i),i=1,nl) |
|---|
| 419 | c sb3d print *,'ICB ',icb |
|---|
| 420 | C |
|---|
| 421 | C *** MAKE SURE THAT COLUMN IS DRY ADIABATIC BETWEEN THE SURFACE *** |
|---|
| 422 | C *** AND CLOUD BASE, AND THAT LIFTED AIR IS POSITIVELY BUOYANT *** |
|---|
| 423 | C *** AT CLOUD BASE *** |
|---|
| 424 | C *** IF NOT, RETURN TO CALLING PROGRAM AFTER RESETTING *** |
|---|
| 425 | C *** SIG(I) AND W0(I) *** |
|---|
| 426 | C |
|---|
| 427 | Cjyg |
|---|
| 428 | CCC TDIF=TVP(ICB)-TV(ICB) |
|---|
| 429 | TDIF = BUOY(ICB) |
|---|
| 430 | ATH1=TH(1) |
|---|
| 431 | Cjyg |
|---|
| 432 | CCC ATH=TH(ICB-1)-1.0 |
|---|
| 433 | ATH=TH(ICB-1)-5.0 |
|---|
| 434 | c ATH=0. ! ajout sb |
|---|
| 435 | c IF (ICB.GT.1) ATH=TH(ICB-1)-5.0 ! modif sb |
|---|
| 436 | IF(TDIF.LT.DTCRIT.OR.ATH.GT.ATH1)THEN |
|---|
| 437 | DO 60 I=1,NL |
|---|
| 438 | SIG(I)=BETA*SIG(I)-2.*ALPHA*TDIF*TDIF |
|---|
| 439 | SIG(I)=AMAX1(SIG(I),0.0) |
|---|
| 440 | W0(I)=BETA*W0(I) |
|---|
| 441 | 60 CONTINUE |
|---|
| 442 | IFLAG=0 |
|---|
| 443 | RETURN |
|---|
| 444 | END IF |
|---|
| 445 | C |
|---|
| 446 | |
|---|
| 447 | |
|---|
| 448 | C *** IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY *** |
|---|
| 449 | C *** NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS *** |
|---|
| 450 | C |
|---|
| 451 | DO 70 I=1,NL |
|---|
| 452 | HP(I)=H(I) |
|---|
| 453 | WT(I)=0.001 |
|---|
| 454 | RP(I)=RR(I) |
|---|
| 455 | UP(I)=U(I) |
|---|
| 456 | VP(I)=V(I) |
|---|
| 457 | DO 71 J=1,NTRA |
|---|
| 458 | TRAP(I,J)=TRA(I,J) |
|---|
| 459 | 71 CONTINUE |
|---|
| 460 | NENT(I)=0 |
|---|
| 461 | WATER(I)=0.0 |
|---|
| 462 | EVAP(I)=0.0 |
|---|
| 463 | B(I)=0.0 |
|---|
| 464 | MP(I)=0.0 |
|---|
| 465 | M(I)=0.0 |
|---|
| 466 | LV(I)=ALV0-CPVMCL*(T(I)-273.15) |
|---|
| 467 | LVCP(I)=LV(I)/CPN(I) |
|---|
| 468 | DO 70 J=1,NL |
|---|
| 469 | QENT(I,J)=RR(J) |
|---|
| 470 | ELIJ(I,J)=0.0 |
|---|
| 471 | MENT(I,J)=0.0 |
|---|
| 472 | SIJ(I,J)=0.0 |
|---|
| 473 | UENT(I,J)=U(J) |
|---|
| 474 | VENT(I,J)=V(J) |
|---|
| 475 | DO 70 K=1,NTRA |
|---|
| 476 | TRAENT(I,J,K)=TRA(J,K) |
|---|
| 477 | 70 CONTINUE |
|---|
| 478 | |
|---|
| 479 | DELTI=1.0/DELT |
|---|
| 480 | C |
|---|
| 481 | C *** FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S *** |
|---|
| 482 | C *** LEVEL OF NEUTRAL BUOYANCY *** |
|---|
| 483 | C |
|---|
| 484 | INB=NL-1 |
|---|
| 485 | DO 80 I=ICB,NL-1 |
|---|
| 486 | Cjyg |
|---|
| 487 | CCC IF((TVP(I)-TV(I)).LT.DTCRIT)THEN |
|---|
| 488 | IF(BUOY(I).LT.DTOVSH)THEN |
|---|
| 489 | INB=MIN(INB,I) |
|---|
| 490 | END IF |
|---|
| 491 | 80 CONTINUE |
|---|
| 492 | |
|---|
| 493 | |
|---|
| 494 | |
|---|
| 495 | |
|---|
| 496 | C |
|---|
| 497 | C *** RESET SIG(I) AND W0(I) FOR I>INB AND I<ICB *** |
|---|
| 498 | C |
|---|
| 499 | IF(INB.LT.(NL-1))THEN |
|---|
| 500 | DO 85 I=INB+1,NL-1 |
|---|
| 501 | Cjyg |
|---|
| 502 | CCC SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(INB)-TVP(INB))* |
|---|
| 503 | CCC 1 ABS(TV(INB)-TVP(INB)) |
|---|
| 504 | SIG(I)=BETA*SIG(I)+2.*ALPHA*BUOY(INB)* |
|---|
| 505 | 1 ABS(BUOY(INB)) |
|---|
| 506 | SIG(I)=AMAX1(SIG(I),0.0) |
|---|
| 507 | W0(I)=BETA*W0(I) |
|---|
| 508 | 85 CONTINUE |
|---|
| 509 | END IF |
|---|
| 510 | DO 87 I=1,ICB |
|---|
| 511 | Cjyg |
|---|
| 512 | CCC SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(ICB)-TVP(ICB))* |
|---|
| 513 | CCC 1 (TV(ICB)-TVP(ICB)) |
|---|
| 514 | SIG(I)=BETA*SIG(I)-2.*ALPHA*BUOY(ICB)*BUOY(ICB) |
|---|
| 515 | SIG(I)=AMAX1(SIG(I),0.0) |
|---|
| 516 | W0(I)=BETA*W0(I) |
|---|
| 517 | 87 CONTINUE |
|---|
| 518 | C |
|---|
| 519 | C *** RESET FRACTIONAL AREAS OF UPDRAFTS AND W0 AT INITIAL TIME *** |
|---|
| 520 | C *** AND AFTER 10 TIME STEPS OF NO CONVECTION *** |
|---|
| 521 | C |
|---|
| 522 | |
|---|
| 523 | IF(SIG(ND).LT.1.5.OR.SIG(ND).GT.12.0)THEN |
|---|
| 524 | DO 90 I=1,NL-1 |
|---|
| 525 | SIG(I)=0.0 |
|---|
| 526 | W0(I)=0.0 |
|---|
| 527 | 90 CONTINUE |
|---|
| 528 | END IF |
|---|
| 529 | C |
|---|
| 530 | C *** CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL *** |
|---|
| 531 | C |
|---|
| 532 | DO 95 I=ICB,INB |
|---|
| 533 | HP(I)=H(1)+(LV(I)+(CPD-CPV)*T(I))*EP(I)*CLW(I) |
|---|
| 534 | 95 CONTINUE |
|---|
| 535 | C |
|---|
| 536 | C *** CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY (CAPE), *** |
|---|
| 537 | C *** VERTICAL VELOCITY (W), FRACTIONAL AREA COVERED BY *** |
|---|
| 538 | C *** UNDILUTE UPDRAFT (SIG), AND UPDRAFT MASS FLUX (M) *** |
|---|
| 539 | C |
|---|
| 540 | CAPE=0.0 |
|---|
| 541 | C |
|---|
| 542 | DO 98 I=ICB+1,INB |
|---|
| 543 | Cjyg1 |
|---|
| 544 | CCC CAPE=CAPE+RD*(TVP(I-1)-TV(I-1))*(PH(I-1)-PH(I))/P(I-1) |
|---|
| 545 | CCC DCAPE=RD*BUOY(I-1)*(PH(I-1)-PH(I))/P(I-1) |
|---|
| 546 | CCC DLNP=(PH(I-1)-PH(I))/P(I-1) |
|---|
| 547 | C The interval on which CAPE is computed starts at PBASE : |
|---|
| 548 | DELTAP = MIN(PBASE,PH(I-1))-MIN(PBASE,PH(I)) |
|---|
| 549 | CAPE=CAPE+RD*BUOY(I-1)*DELTAP/P(I-1) |
|---|
| 550 | DCAPE=RD*BUOY(I-1)*DELTAP/P(I-1) |
|---|
| 551 | DLNP=DELTAP/P(I-1) |
|---|
| 552 | Cjyg2 |
|---|
| 553 | c sb3d print *,'buoy,dlnp,dcape,cape',buoy(i-1),dlnp,dcape,cape |
|---|
| 554 | CAPE=AMAX1(0.0,CAPE) |
|---|
| 555 | C |
|---|
| 556 | SIGOLD=SIG(I) |
|---|
| 557 | DTMIN=100.0 |
|---|
| 558 | DO 97 J=ICB,I-1 |
|---|
| 559 | Cjyg |
|---|
| 560 | CCC DTMIN=AMIN1(DTMIN,(TVP(J)-TV(J))) |
|---|
| 561 | DTMIN=AMIN1(DTMIN,BUOY(J)) |
|---|
| 562 | 97 CONTINUE |
|---|
| 563 | c sb3d print *, 'DTMIN, BETA, ALPHA, SIG = ',DTMIN,BETA,ALPHA,SIG(I) |
|---|
| 564 | SIG(I)=BETA*SIG(I)+ALPHA*DTMIN*ABS(DTMIN) |
|---|
| 565 | SIG(I)=AMAX1(SIG(I),0.0) |
|---|
| 566 | SIG(I)=AMIN1(SIG(I),0.01) |
|---|
| 567 | FAC=AMIN1(((DTCRIT-DTMIN)/DTCRIT),1.0) |
|---|
| 568 | Cjyg |
|---|
| 569 | CC Essais de reduction de la vitesse |
|---|
| 570 | CC FAC = FAC*.5 |
|---|
| 571 | C |
|---|
| 572 | W=(1.-BETA)*FAC*SQRT(CAPE)+BETA*W0(I) |
|---|
| 573 | AMU=0.5*(SIG(I)+SIGOLD)*W |
|---|
| 574 | M(I)=AMU*0.007*P(I)*(PH(I)-PH(I+1))/TV(I) |
|---|
| 575 | W0(I)=W |
|---|
| 576 | 98 CONTINUE |
|---|
| 577 | W0(ICB)=0.5*W0(ICB+1) |
|---|
| 578 | M(ICB)=0.5*M(ICB+1)*(PH(ICB)-PH(ICB+1))/(PH(ICB+1)-PH(ICB+2)) |
|---|
| 579 | SIG(ICB)=SIG(ICB+1) |
|---|
| 580 | SIG(ICB-1)=SIG(ICB) |
|---|
| 581 | cjyg1 |
|---|
| 582 | c sb3d print *, 'Cloud base, c. top, CAPE',ICB,INB,cape |
|---|
| 583 | c sb3d print *, 'SIG ',(sig(i),i=1,inb) |
|---|
| 584 | c sb3d print *, 'W ',(w0(i),i=1,inb) |
|---|
| 585 | c sb3d print *, 'M ',(m(i), i=1,inb) |
|---|
| 586 | c sb3d print *, 'Dt1 ',(tvp(i)-tv(i),i=1,inb) |
|---|
| 587 | c sb3d print *, 'Dt_vrai ',(buoy(i),i=1,inb) |
|---|
| 588 | Cjyg2 |
|---|
| 589 | C |
|---|
| 590 | C *** CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING *** |
|---|
| 591 | C *** RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING *** |
|---|
| 592 | C *** FRACTION (SIJ) *** |
|---|
| 593 | C |
|---|
| 594 | |
|---|
| 595 | |
|---|
| 596 | DO 170 I=ICB,INB |
|---|
| 597 | RTI=RR(1)-EP(I)*CLW(I) |
|---|
| 598 | DO 160 J=ICB-1,INB |
|---|
| 599 | BF2=1.+LV(J)*LV(J)*RS(J)/(RV*T(J)*T(J)*CPD) |
|---|
| 600 | ANUM=H(J)-HP(I)+(CPV-CPD)*T(J)*(RTI-RR(J)) |
|---|
| 601 | DENOM=H(I)-HP(I)+(CPD-CPV)*(RR(I)-RTI)*T(J) |
|---|
| 602 | DEI=DENOM |
|---|
| 603 | IF(ABS(DEI).LT.0.01)DEI=0.01 |
|---|
| 604 | SIJ(I,J)=ANUM/DEI |
|---|
| 605 | SIJ(I,I)=1.0 |
|---|
| 606 | ALTEM=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI-RS(J) |
|---|
| 607 | ALTEM=ALTEM/BF2 |
|---|
| 608 | CWAT=CLW(J)*(1.-EP(J)) |
|---|
| 609 | STEMP=SIJ(I,J) |
|---|
| 610 | IF((STEMP.LT.0.0.OR.STEMP.GT.1.0.OR. |
|---|
| 611 | 1 ALTEM.GT.CWAT).AND.J.GT.I)THEN |
|---|
| 612 | ANUM=ANUM-LV(J)*(RTI-RS(J)-CWAT*BF2) |
|---|
| 613 | DENOM=DENOM+LV(J)*(RR(I)-RTI) |
|---|
| 614 | IF(ABS(DENOM).LT.0.01)DENOM=0.01 |
|---|
| 615 | SIJ(I,J)=ANUM/DENOM |
|---|
| 616 | ALTEM=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI-RS(J) |
|---|
| 617 | ALTEM=ALTEM-(BF2-1.)*CWAT |
|---|
| 618 | END IF |
|---|
| 619 | |
|---|
| 620 | |
|---|
| 621 | IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.95)THEN |
|---|
| 622 | QENT(I,J)=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI |
|---|
| 623 | UENT(I,J)=SIJ(I,J)*U(I)+(1.-SIJ(I,J))*U(NK) |
|---|
| 624 | VENT(I,J)=SIJ(I,J)*V(I)+(1.-SIJ(I,J))*V(NK) |
|---|
| 625 | DO K=1,NTRA |
|---|
| 626 | TRAENT(I,J,K)=SIJ(I,J)*TRA(I,K)+(1.-SIJ(I,J))* |
|---|
| 627 | 1 TRA(NK,K) |
|---|
| 628 | END DO |
|---|
| 629 | ELIJ(I,J)=ALTEM |
|---|
| 630 | ELIJ(I,J)=AMAX1(0.0,ELIJ(I,J)) |
|---|
| 631 | MENT(I,J)=M(I)/(1.-SIJ(I,J)) |
|---|
| 632 | NENT(I)=NENT(I)+1 |
|---|
| 633 | END IF |
|---|
| 634 | SIJ(I,J)=AMAX1(0.0,SIJ(I,J)) |
|---|
| 635 | SIJ(I,J)=AMIN1(1.0,SIJ(I,J)) |
|---|
| 636 | 160 CONTINUE |
|---|
| 637 | C |
|---|
| 638 | C *** IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS *** |
|---|
| 639 | C *** AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES *** |
|---|
| 640 | C |
|---|
| 641 | IF(NENT(I).EQ.0)THEN |
|---|
| 642 | MENT(I,I)=M(I) |
|---|
| 643 | QENT(I,I)=RR(NK)-EP(I)*CLW(I) |
|---|
| 644 | UENT(I,I)=U(NK) |
|---|
| 645 | VENT(I,I)=V(NK) |
|---|
| 646 | DO J=1,NTRA |
|---|
| 647 | TRAENT(I,I,J)=TRA(NK,J) |
|---|
| 648 | END DO |
|---|
| 649 | ELIJ(I,I)=CLW(I) |
|---|
| 650 | SIJ(I,I)=1.0 |
|---|
| 651 | END IF |
|---|
| 652 | C |
|---|
| 653 | DO J = ICB-1,INB |
|---|
| 654 | SIGIJ(I,J) = SIJ(I,J) |
|---|
| 655 | ENDDO |
|---|
| 656 | C |
|---|
| 657 | 170 CONTINUE |
|---|
| 658 | C |
|---|
| 659 | C *** NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL *** |
|---|
| 660 | C *** PROBABILITIES OF MIXING *** |
|---|
| 661 | C |
|---|
| 662 | |
|---|
| 663 | DO 200 I=ICB,INB |
|---|
| 664 | IF(NENT(I).NE.0)THEN |
|---|
| 665 | QP=RR(1)-EP(I)*CLW(I) |
|---|
| 666 | ANUM=H(I)-HP(I)-LV(I)*(QP-RS(I))+(CPV-CPD)*T(I)* |
|---|
| 667 | 1 (QP-RR(I)) |
|---|
| 668 | DENOM=H(I)-HP(I)+LV(I)*(RR(I)-QP)+ |
|---|
| 669 | 1 (CPD-CPV)*T(I)*(RR(I)-QP) |
|---|
| 670 | IF(ABS(DENOM).LT.0.01)DENOM=0.01 |
|---|
| 671 | SCRIT=ANUM/DENOM |
|---|
| 672 | ALT=QP-RS(I)+SCRIT*(RR(I)-QP) |
|---|
| 673 | IF(SCRIT.LE.0.0.OR.ALT.LE.0.0)SCRIT=1.0 |
|---|
| 674 | SMAX=0.0 |
|---|
| 675 | ASIJ=0.0 |
|---|
| 676 | DO 175 J=INB,ICB-1,-1 |
|---|
| 677 | IF(SIJ(I,J).GT.1.0E-16.AND.SIJ(I,J).LT.0.95)THEN |
|---|
| 678 | WGH=1.0 |
|---|
| 679 | IF(J.GT.I)THEN |
|---|
| 680 | SJMAX=AMAX1(SIJ(I,J+1),SMAX) |
|---|
| 681 | SJMAX=AMIN1(SJMAX,SCRIT) |
|---|
| 682 | SMAX=AMAX1(SIJ(I,J),SMAX) |
|---|
| 683 | SJMIN=AMAX1(SIJ(I,J-1),SMAX) |
|---|
| 684 | SJMIN=AMIN1(SJMIN,SCRIT) |
|---|
| 685 | IF(SIJ(I,J).LT.(SMAX-1.0E-16))WGH=0.0 |
|---|
| 686 | SMID=AMIN1(SIJ(I,J),SCRIT) |
|---|
| 687 | ELSE |
|---|
| 688 | SJMAX=AMAX1(SIJ(I,J+1),SCRIT) |
|---|
| 689 | SMID=AMAX1(SIJ(I,J),SCRIT) |
|---|
| 690 | SJMIN=0.0 |
|---|
| 691 | IF(J.GT.1)SJMIN=SIJ(I,J-1) |
|---|
| 692 | SJMIN=AMAX1(SJMIN,SCRIT) |
|---|
| 693 | END IF |
|---|
| 694 | DELP=ABS(SJMAX-SMID) |
|---|
| 695 | DELM=ABS(SJMIN-SMID) |
|---|
| 696 | ASIJ=ASIJ+WGH*(DELP+DELM) |
|---|
| 697 | MENT(I,J)=MENT(I,J)*(DELP+DELM)*WGH |
|---|
| 698 | END IF |
|---|
| 699 | 175 CONTINUE |
|---|
| 700 | ASIJ=AMAX1(1.0E-16,ASIJ) |
|---|
| 701 | ASIJ=1.0/ASIJ |
|---|
| 702 | DO 180 J=ICB-1,INB |
|---|
| 703 | MENT(I,J)=MENT(I,J)*ASIJ |
|---|
| 704 | 180 CONTINUE |
|---|
| 705 | ASUM=0.0 |
|---|
| 706 | BSUM=0.0 |
|---|
| 707 | DO 190 J=ICB-1,INB |
|---|
| 708 | ASUM=ASUM+MENT(I,J) |
|---|
| 709 | MENT(I,J)=MENT(I,J)*SIG(J) |
|---|
| 710 | BSUM=BSUM+MENT(I,J) |
|---|
| 711 | 190 CONTINUE |
|---|
| 712 | BSUM=AMAX1(BSUM,1.0E-16) |
|---|
| 713 | BSUM=1.0/BSUM |
|---|
| 714 | DO 195 J=ICB-1,INB |
|---|
| 715 | MENT(I,J)=MENT(I,J)*ASUM*BSUM |
|---|
| 716 | 195 CONTINUE |
|---|
| 717 | CSUM=0.0 |
|---|
| 718 | DO 197 J=ICB-1,INB |
|---|
| 719 | CSUM=CSUM+MENT(I,J) |
|---|
| 720 | 197 CONTINUE |
|---|
| 721 | |
|---|
| 722 | IF(CSUM.LT.M(I))THEN |
|---|
| 723 | NENT(I)=0 |
|---|
| 724 | MENT(I,I)=M(I) |
|---|
| 725 | QENT(I,I)=RR(1)-EP(I)*CLW(I) |
|---|
| 726 | UENT(I,I)=U(NK) |
|---|
| 727 | VENT(I,I)=V(NK) |
|---|
| 728 | DO J=1,NTRA |
|---|
| 729 | TRAENT(I,I,J)=TRA(NK,J) |
|---|
| 730 | END DO |
|---|
| 731 | ELIJ(I,I)=CLW(I) |
|---|
| 732 | SIJ(I,I)=1.0 |
|---|
| 733 | END IF |
|---|
| 734 | END IF |
|---|
| 735 | 200 CONTINUE |
|---|
| 736 | |
|---|
| 737 | |
|---|
| 738 | |
|---|
| 739 | *************************************************************** |
|---|
| 740 | ** CALCUL DES MENTS(I,J) ET DES QENTS(I,J) |
|---|
| 741 | ************************************************************** |
|---|
| 742 | |
|---|
| 743 | DO im=1,nd |
|---|
| 744 | do jm=1,nd |
|---|
| 745 | |
|---|
| 746 | QENTS(im,jm)=QENT(im,jm) |
|---|
| 747 | MENTS(im,jm)=MENT(im,jm) |
|---|
| 748 | enddo |
|---|
| 749 | enddo |
|---|
| 750 | |
|---|
| 751 | *********************************************************** |
|---|
| 752 | |
|---|
| 753 | |
|---|
| 754 | |
|---|
| 755 | C |
|---|
| 756 | C *** CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING *** |
|---|
| 757 | C *** DOWNDRAFT CALCULATION *** |
|---|
| 758 | C |
|---|
| 759 | IF(EP(INB).LT.0.0001)GOTO 405 |
|---|
| 760 | C |
|---|
| 761 | C *** INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER *** |
|---|
| 762 | C *** AND CONDENSED WATER FLUX *** |
|---|
| 763 | C |
|---|
| 764 | WFLUX=0.0 |
|---|
| 765 | TINV=1./3. |
|---|
| 766 | C |
|---|
| 767 | C *** BEGIN DOWNDRAFT LOOP *** |
|---|
| 768 | C |
|---|
| 769 | DO 400 I=INB,1,-1 |
|---|
| 770 | C |
|---|
| 771 | C *** CALCULATE DETRAINED PRECIPITATION *** |
|---|
| 772 | C |
|---|
| 773 | |
|---|
| 774 | |
|---|
| 775 | WDTRAIN=10.0*EP(I)*M(I)*CLW(I) |
|---|
| 776 | IF(I.GT.1)THEN |
|---|
| 777 | DO 320 J=1,I-1 |
|---|
| 778 | AWAT=ELIJ(J,I)-(1.-EP(I))*CLW(I) |
|---|
| 779 | AWAT=AMAX1(AWAT,0.0) |
|---|
| 780 | 320 WDTRAIN=WDTRAIN+10.0*AWAT*MENT(J,I) |
|---|
| 781 | END IF |
|---|
| 782 | C |
|---|
| 783 | C *** FIND RAIN WATER AND EVAPORATION USING PROVISIONAL *** |
|---|
| 784 | C *** ESTIMATES OF RP(I)AND RP(I-1) *** |
|---|
| 785 | C |
|---|
| 786 | |
|---|
| 787 | |
|---|
| 788 | WT(I)=45.0 |
|---|
| 789 | IF(I.LT.INB)THEN |
|---|
| 790 | RP(I)=RP(I+1)+(CPD*(T(I+1)-T(I))+GZ(I+1)-GZ(I))/LV(I) |
|---|
| 791 | RP(I)=0.5*(RP(I)+RR(I)) |
|---|
| 792 | END IF |
|---|
| 793 | RP(I)=AMAX1(RP(I),0.0) |
|---|
| 794 | RP(I)=AMIN1(RP(I),RS(I)) |
|---|
| 795 | RP(INB)=RR(INB) |
|---|
| 796 | IF(I.EQ.1)THEN |
|---|
| 797 | AFAC=P(1)*(RS(1)-RP(1))/(1.0E4+2000.0*P(1)*RS(1)) |
|---|
| 798 | ELSE |
|---|
| 799 | RP(I-1)=RP(I)+(CPD*(T(I)-T(I-1))+GZ(I)-GZ(I-1))/LV(I) |
|---|
| 800 | RP(I-1)=0.5*(RP(I-1)+RR(I-1)) |
|---|
| 801 | RP(I-1)=AMIN1(RP(I-1),RS(I-1)) |
|---|
| 802 | RP(I-1)=AMAX1(RP(I-1),0.0) |
|---|
| 803 | AFAC1=P(I)*(RS(I)-RP(I))/(1.0E4+2000.0*P(I)*RS(I)) |
|---|
| 804 | AFAC2=P(I-1)*(RS(I-1)-RP(I-1))/(1.0E4+ |
|---|
| 805 | 1 2000.0*P(I-1)*RS(I-1)) |
|---|
| 806 | AFAC=0.5*(AFAC1+AFAC2) |
|---|
| 807 | END IF |
|---|
| 808 | IF(I.EQ.INB)AFAC=0.0 |
|---|
| 809 | AFAC=AMAX1(AFAC,0.0) |
|---|
| 810 | BFAC=1./(SIGD*WT(I)) |
|---|
| 811 | C |
|---|
| 812 | Cjyg1 |
|---|
| 813 | CCC SIGT=1.0 |
|---|
| 814 | CCC IF(I.GE.ICB)SIGT=SIGP(I) |
|---|
| 815 | C Prise en compte de la variation progressive de SIGT dans |
|---|
| 816 | C les couches ICB et ICB-1: |
|---|
| 817 | C Pour PLCL<PH(I+1), PR1=0 & PR2=1 |
|---|
| 818 | C Pour PLCL>PH(I), PR1=1 & PR2=0 |
|---|
| 819 | C Pour PH(I+1)<PLCL<PH(I), PR1 est la proportion a cheval |
|---|
| 820 | C sur le nuage, et PR2 est la proportion sous la base du |
|---|
| 821 | C nuage. |
|---|
| 822 | PR1 =(PLCL-PH(I+1))/(PH(I)-PH(I+1)) |
|---|
| 823 | PR1 = MAX(0.,MIN(1.,PR1)) |
|---|
| 824 | PR2 = (PH(I)-PLCL)/(PH(I)-PH(I+1)) |
|---|
| 825 | PR2 = MAX(0.,MIN(1.,PR2)) |
|---|
| 826 | SIGT = SIGP(I)*PR1 + PR2 |
|---|
| 827 | c sb3d print *,'i,sigt,pr1,pr2', i,sigt,pr1,pr2 |
|---|
| 828 | Cjyg2 |
|---|
| 829 | C |
|---|
| 830 | B6=BFAC*50.*SIGD*(PH(I)-PH(I+1))*SIGT*AFAC |
|---|
| 831 | C6=WATER(I+1)+BFAC*WDTRAIN-50.*SIGD*BFAC* |
|---|
| 832 | 1 (PH(I)-PH(I+1))*EVAP(I+1) |
|---|
| 833 | IF(C6.GT.0.0)THEN |
|---|
| 834 | REVAP=0.5*(-B6+SQRT(B6*B6+4.*C6)) |
|---|
| 835 | EVAP(I)=SIGT*AFAC*REVAP |
|---|
| 836 | WATER(I)=REVAP*REVAP |
|---|
| 837 | ELSE |
|---|
| 838 | EVAP(I)=-EVAP(I+1)+0.02*(WDTRAIN+SIGD*WT(I)* |
|---|
| 839 | 1 WATER(I+1))/(SIGD*(PH(I)-PH(I+1))) |
|---|
| 840 | END IF |
|---|
| 841 | |
|---|
| 842 | |
|---|
| 843 | C |
|---|
| 844 | C *** CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER *** |
|---|
| 845 | C *** HYDROSTATIC APPROXIMATION *** |
|---|
| 846 | C |
|---|
| 847 | IF(I.EQ.1)GOTO 360 |
|---|
| 848 | TEVAP=AMAX1(0.0,EVAP(I)) |
|---|
| 849 | DELTH=AMAX1(0.001,(TH(I)-TH(I-1))) |
|---|
| 850 | MP(I)=10.*LVCP(I)*SIGD*TEVAP*(P(I-1)-P(I))/DELTH |
|---|
| 851 | C |
|---|
| 852 | C *** IF HYDROSTATIC ASSUMPTION FAILS, *** |
|---|
| 853 | C *** SOLVE CUBIC DIFFERENCE EQUATION FOR DOWNDRAFT THETA *** |
|---|
| 854 | C *** AND MASS FLUX FROM TWO SIMULTANEOUS DIFFERENTIAL EQNS *** |
|---|
| 855 | C |
|---|
| 856 | AMFAC=SIGD*SIGD*70.0*PH(I)*(P(I-1)-P(I))* |
|---|
| 857 | 1 (TH(I)-TH(I-1))/(TV(I)*TH(I)) |
|---|
| 858 | AMP2=ABS(MP(I+1)*MP(I+1)-MP(I)*MP(I)) |
|---|
| 859 | IF(AMP2.GT.(0.1*AMFAC))THEN |
|---|
| 860 | XF=100.0*SIGD*SIGD*SIGD*(PH(I)-PH(I+1)) |
|---|
| 861 | TF=B(I)-5.0*(TH(I)-TH(I-1))*T(I)/(LVCP(I)*SIGD*TH(I)) |
|---|
| 862 | AF=XF*TF+MP(I+1)*MP(I+1)*TINV |
|---|
| 863 | BF=2.*(TINV*MP(I+1))**3+TINV*MP(I+1)*XF*TF+50.* |
|---|
| 864 | 1 (P(I-1)-P(I))*XF*TEVAP |
|---|
| 865 | FAC2=1.0 |
|---|
| 866 | IF(BF.LT.0.0)FAC2=-1.0 |
|---|
| 867 | BF=ABS(BF) |
|---|
| 868 | UR=0.25*BF*BF-AF*AF*AF*TINV*TINV*TINV |
|---|
| 869 | IF(UR.GE.0.0)THEN |
|---|
| 870 | SRU=SQRT(UR) |
|---|
| 871 | FAC=1.0 |
|---|
| 872 | IF((0.5*BF-SRU).LT.0.0)FAC=-1.0 |
|---|
| 873 | MP(I)=MP(I+1)*TINV+(0.5*BF+SRU)**TINV+ |
|---|
| 874 | 1 FAC*(ABS(0.5*BF-SRU))**TINV |
|---|
| 875 | ELSE |
|---|
| 876 | D=ATAN(2.*SQRT(-UR)/(BF+1.0E-28)) |
|---|
| 877 | IF(FAC2.LT.0.0)D=3.14159-D |
|---|
| 878 | MP(I)=MP(I+1)*TINV+2.*SQRT(AF*TINV)*COS(D*TINV) |
|---|
| 879 | END IF |
|---|
| 880 | MP(I)=AMAX1(0.0,MP(I)) |
|---|
| 881 | B(I-1)=B(I)+100.0*(P(I-1)-P(I))*TEVAP/(MP(I)+SIGD*0.1)- |
|---|
| 882 | 1 10.0*(TH(I)-TH(I-1))*T(I)/(LVCP(I)*SIGD*TH(I)) |
|---|
| 883 | B(I-1)=AMAX1(B(I-1),0.0) |
|---|
| 884 | END IF |
|---|
| 885 | |
|---|
| 886 | |
|---|
| 887 | C |
|---|
| 888 | C *** LIMIT MAGNITUDE OF MP(I) TO MEET CFL CONDITION *** |
|---|
| 889 | C |
|---|
| 890 | AMPMAX=2.0*(PH(I)-PH(I+1))*DELTI |
|---|
| 891 | AMP2=2.0*(PH(I-1)-PH(I))*DELTI |
|---|
| 892 | AMPMAX=AMIN1(AMPMAX,AMP2) |
|---|
| 893 | MP(I)=AMIN1(MP(I),AMPMAX) |
|---|
| 894 | C |
|---|
| 895 | C *** FORCE MP TO DECREASE LINEARLY TO ZERO *** |
|---|
| 896 | C *** BETWEEN CLOUD BASE AND THE SURFACE *** |
|---|
| 897 | C |
|---|
| 898 | IF(P(I).GT.P(ICB))THEN |
|---|
| 899 | MP(I)=MP(ICB)*(P(1)-P(I))/(P(1)-P(ICB)) |
|---|
| 900 | END IF |
|---|
| 901 | 360 CONTINUE |
|---|
| 902 | C |
|---|
| 903 | C *** FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT *** |
|---|
| 904 | C |
|---|
| 905 | IF(I.EQ.INB)GOTO 400 |
|---|
| 906 | RP(I)=RR(I) |
|---|
| 907 | IF(MP(I).GT.MP(I+1))THEN |
|---|
| 908 | RP(I)=RP(I+1)*MP(I+1)+RR(I)*(MP(I)-MP(I+1))+ |
|---|
| 909 | 1 5.*SIGD*(PH(I)-PH(I+1))*(EVAP(I+1)+EVAP(I)) |
|---|
| 910 | RP(I)=RP(I)/MP(I) |
|---|
| 911 | UP(I)=UP(I+1)*MP(I+1)+U(I)*(MP(I)-MP(I+1)) |
|---|
| 912 | UP(I)=UP(I)/MP(I) |
|---|
| 913 | VP(I)=VP(I+1)*MP(I+1)+V(I)*(MP(I)-MP(I+1)) |
|---|
| 914 | VP(I)=VP(I)/MP(I) |
|---|
| 915 | DO J=1,NTRA |
|---|
| 916 | TRAP(I,J)=TRAP(I+1,J)*MP(I+1)+ |
|---|
| 917 | s TRAP(I,J)*(MP(I)-MP(I+1)) |
|---|
| 918 | TRAP(I,J)=TRAP(I,J)/MP(I) |
|---|
| 919 | END DO |
|---|
| 920 | ELSE |
|---|
| 921 | IF(MP(I+1).GT.1.0E-16)THEN |
|---|
| 922 | RP(I)=RP(I+1)+5.0*SIGD*(PH(I)-PH(I+1))*(EVAP(I+1)+ |
|---|
| 923 | 1 EVAP(I))/MP(I+1) |
|---|
| 924 | UP(I)=UP(I+1) |
|---|
| 925 | VP(I)=VP(I+1) |
|---|
| 926 | DO J=1,NTRA |
|---|
| 927 | TRAP(I,J)=TRAP(I+1,J) |
|---|
| 928 | END DO |
|---|
| 929 | END IF |
|---|
| 930 | END IF |
|---|
| 931 | RP(I)=AMIN1(RP(I),RS(I)) |
|---|
| 932 | RP(I)=AMAX1(RP(I),0.0) |
|---|
| 933 | 400 CONTINUE |
|---|
| 934 | C |
|---|
| 935 | C *** CALCULATE SURFACE PRECIPITATION IN MM/DAY *** |
|---|
| 936 | C |
|---|
| 937 | PRECIP=WT(1)*SIGD*WATER(1)*8640.0 |
|---|
| 938 | 405 CONTINUE |
|---|
| 939 | C |
|---|
| 940 | C *** CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE *** |
|---|
| 941 | C *** AND MIXING RATIO *** |
|---|
| 942 | C |
|---|
| 943 | DPINV=1.0/(PH(1)-PH(2)) |
|---|
| 944 | AM=0.0 |
|---|
| 945 | DO 410 K=2,INB |
|---|
| 946 | 410 AM=AM+M(K) |
|---|
| 947 | IF((0.1*DPINV*AM).GE.DELTI)IFLAG=4 |
|---|
| 948 | FT(1)=0.1*DPINV*AM*(T(2)-T(1)+(GZ(2)-GZ(1))/CPN(1)) |
|---|
| 949 | FT(1)=FT(1)-0.5*LVCP(1)*SIGD*(EVAP(1)+EVAP(2)) |
|---|
| 950 | FT(1)=FT(1)-0.09*SIGD*MP(2)*T(1)*B(1)*DPINV |
|---|
| 951 | FT(1)=FT(1)+0.01*SIGD*WT(1)*(CL-CPD)*WATER(2)*(T(2)- |
|---|
| 952 | 1 T(1))*DPINV/CPN(1) |
|---|
| 953 | FR(1)=0.1*MP(2)*(RP(2)-RR(1))* |
|---|
| 954 | 1 DPINV+SIGD*0.5*(EVAP(1)+EVAP(2)) |
|---|
| 955 | FR(1)=FR(1)+0.1*AM*(RR(2)-RR(1))*DPINV |
|---|
| 956 | FU(1)=FU(1)+0.1*DPINV*(MP(2)*(UP(2)-U(1))+AM*(U(2)-U(1))) |
|---|
| 957 | FV(1)=FV(1)+0.1*DPINV*(MP(2)*(VP(2)-V(1))+AM*(V(2)-V(1))) |
|---|
| 958 | DO J=1,NTRA |
|---|
| 959 | FTRA(1,J)=FTRA(1,J)+0.1*DPINV*(MP(2)*(TRAP(2,J)-TRA(1,J))+ |
|---|
| 960 | 1 AM*(TRA(2,J)-TRA(1,J))) |
|---|
| 961 | END DO |
|---|
| 962 | AMDE=0.0 |
|---|
| 963 | DO 415 J=2,INB |
|---|
| 964 | FR(1)=FR(1)+0.1*DPINV*MENT(J,1)*(QENT(J,1)-RR(1)) |
|---|
| 965 | FU(1)=FU(1)+0.1*DPINV*MENT(J,1)*(UENT(J,1)-U(1)) |
|---|
| 966 | FV(1)=FV(1)+0.1*DPINV*MENT(J,1)*(VENT(J,1)-V(1)) |
|---|
| 967 | DO K=1,NTRA |
|---|
| 968 | FTRA(1,K)=FTRA(1,K)+0.1*DPINV*MENT(J,1)*(TRAENT(J,1,K)- |
|---|
| 969 | 1 TRA(1,K)) |
|---|
| 970 | END DO |
|---|
| 971 | 415 CONTINUE |
|---|
| 972 | C |
|---|
| 973 | C *** CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO *** |
|---|
| 974 | C *** AT LEVELS ABOVE THE LOWEST LEVEL *** |
|---|
| 975 | C |
|---|
| 976 | C *** FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES *** |
|---|
| 977 | C *** THROUGH EACH LEVEL *** |
|---|
| 978 | C |
|---|
| 979 | |
|---|
| 980 | |
|---|
| 981 | DO 500 I=2,INB |
|---|
| 982 | DPINV=1.0/(PH(I)-PH(I+1)) |
|---|
| 983 | CPINV=1.0/CPN(I) |
|---|
| 984 | AMP1=0.0 |
|---|
| 985 | DO 440 K=I+1,INB+1 |
|---|
| 986 | 440 AMP1=AMP1+M(K) |
|---|
| 987 | DO 450 K=1,I |
|---|
| 988 | DO 450 J=I+1,INB+1 |
|---|
| 989 | AMP1=AMP1+MENT(K,J) |
|---|
| 990 | 450 CONTINUE |
|---|
| 991 | IF((0.1*DPINV*AMP1).GE.DELTI)IFLAG=4 |
|---|
| 992 | AD=0.0 |
|---|
| 993 | DO 470 K=1,I-1 |
|---|
| 994 | DO 470 J=I,INB |
|---|
| 995 | 470 AD=AD+MENT(J,K) |
|---|
| 996 | FT(I)=0.1*DPINV*(AMP1*(T(I+1)-T(I)+(GZ(I+1)-GZ(I))* |
|---|
| 997 | 1 CPINV)-AD*(T(I)-T(I-1)+(GZ(I)-GZ(I-1))*CPINV)) |
|---|
| 998 | 2 -0.5*SIGD*LVCP(I)*(EVAP(I)+EVAP(I+1)) |
|---|
| 999 | RAT=CPN(I-1)*CPINV |
|---|
| 1000 | FT(I)=FT(I)-0.09*SIGD*(MP(I+1)*T(I)* |
|---|
| 1001 | 1 B(I)-MP(I)*T(I-1)*RAT*B(I-1))*DPINV |
|---|
| 1002 | FT(I)=FT(I)+0.1*DPINV*MENT(I,I)*(HP(I)-H(I)+ |
|---|
| 1003 | 1 T(I)*(CPV-CPD)*(RR(I)-QENT(I,I)))*CPINV |
|---|
| 1004 | FT(I)=FT(I)+0.01*SIGD*WT(I)*(CL-CPD)*WATER(I+1)* |
|---|
| 1005 | 1 (T(I+1)-T(I))*DPINV*CPINV |
|---|
| 1006 | FR(I)=0.1*DPINV*(AMP1*(RR(I+1)-RR(I))- |
|---|
| 1007 | 1 AD*(RR(I)-RR(I-1))) |
|---|
| 1008 | FU(I)=FU(I)+0.1*DPINV*(AMP1*(U(I+1)-U(I))- |
|---|
| 1009 | 1 AD*(U(I)-U(I-1))) |
|---|
| 1010 | FV(I)=FV(I)+0.1*DPINV*(AMP1*(V(I+1)-V(I))- |
|---|
| 1011 | 1 AD*(V(I)-V(I-1))) |
|---|
| 1012 | DO K=1,NTRA |
|---|
| 1013 | FTRA(I,K)=FTRA(I,K)+0.1*DPINV*(AMP1*(TRA(I+1,K)- |
|---|
| 1014 | 1 TRA(I,K))-AD*(TRA(I,K)-TRA(I-1,K))) |
|---|
| 1015 | END DO |
|---|
| 1016 | DO 480 K=1,I-1 |
|---|
| 1017 | AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I) |
|---|
| 1018 | AWAT=AMAX1(AWAT,0.0) |
|---|
| 1019 | FR(I)=FR(I)+0.1*DPINV*MENT(K,I)*(QENT(K,I)-AWAT |
|---|
| 1020 | 1 -RR(I)) |
|---|
| 1021 | FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I)) |
|---|
| 1022 | FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I)) |
|---|
| 1023 | DO J=1,NTRA |
|---|
| 1024 | FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)- |
|---|
| 1025 | 1 TRA(I,J)) |
|---|
| 1026 | END DO |
|---|
| 1027 | 480 CONTINUE |
|---|
| 1028 | DO 490 K=I,INB |
|---|
| 1029 | FR(I)=FR(I)+0.1*DPINV*MENT(K,I)*(QENT(K,I)-RR(I)) |
|---|
| 1030 | FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I)) |
|---|
| 1031 | FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I)) |
|---|
| 1032 | DO J=1,NTRA |
|---|
| 1033 | FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)- |
|---|
| 1034 | 1 TRA(I,J)) |
|---|
| 1035 | END DO |
|---|
| 1036 | 490 CONTINUE |
|---|
| 1037 | FR(I)=FR(I)+0.5*SIGD*(EVAP(I)+EVAP(I+1))+0.1*(MP(I+1)* |
|---|
| 1038 | 1 (RP(I+1)-RR(I))-MP(I)*(RP(I)-RR(I-1)))*DPINV |
|---|
| 1039 | FU(I)=FU(I)+0.1*(MP(I+1)*(UP(I+1)-U(I))-MP(I)* |
|---|
| 1040 | 1 (UP(I)-U(I-1)))*DPINV |
|---|
| 1041 | FV(I)=FV(I)+0.1*(MP(I+1)*(VP(I+1)-V(I))-MP(I)* |
|---|
| 1042 | 1 (VP(I)-V(I-1)))*DPINV |
|---|
| 1043 | DO J=1,NTRA |
|---|
| 1044 | FTRA(I,J)=FTRA(I,J)+0.1*DPINV*(MP(I+1)*(TRAP(I+1,J)-TRA(I,J))- |
|---|
| 1045 | 1 MP(I)*(TRAP(I,J)-TRAP(I-1,J))) |
|---|
| 1046 | END DO |
|---|
| 1047 | 500 CONTINUE |
|---|
| 1048 | |
|---|
| 1049 | |
|---|
| 1050 | |
|---|
| 1051 | C |
|---|
| 1052 | C *** MOVE THE DETRAINMENT AT LEVEL INB DOWN TO LEVEL INB-1 *** |
|---|
| 1053 | C *** IN SUCH A WAY AS TO PRESERVE THE VERTICALLY *** |
|---|
| 1054 | C *** INTEGRATED ENTHALPY AND WATER TENDENCIES *** |
|---|
| 1055 | C |
|---|
| 1056 | AX=0.1*MENT(INB,INB)*(HP(INB)-H(INB)+T(INB)* |
|---|
| 1057 | 1 (CPV-CPD)*(RR(INB)-QENT(INB,INB)))/(CPN(INB)* |
|---|
| 1058 | 2 (PH(INB)-PH(INB+1))) |
|---|
| 1059 | FT(INB)=FT(INB)-AX |
|---|
| 1060 | FT(INB-1)=FT(INB-1)+AX*CPN(INB)*(PH(INB)-PH(INB+1))/ |
|---|
| 1061 | 1 (CPN(INB-1)*(PH(INB-1)-PH(INB))) |
|---|
| 1062 | BX=0.1*MENT(INB,INB)*(QENT(INB,INB)-RR(INB))/ |
|---|
| 1063 | 1 (PH(INB)-PH(INB+1)) |
|---|
| 1064 | FR(INB)=FR(INB)-BX |
|---|
| 1065 | FR(INB-1)=FR(INB-1)+BX*(PH(INB)-PH(INB+1))/ |
|---|
| 1066 | 1 (PH(INB-1)-PH(INB)) |
|---|
| 1067 | CX=0.1*MENT(INB,INB)*(UENT(INB,INB)-U(INB))/ |
|---|
| 1068 | 1 (PH(INB)-PH(INB+1)) |
|---|
| 1069 | FU(INB)=FU(INB)-CX |
|---|
| 1070 | FU(INB-1)=FU(INB-1)+CX*(PH(INB)-PH(INB+1))/ |
|---|
| 1071 | 1 (PH(INB-1)-PH(INB)) |
|---|
| 1072 | DX=0.1*MENT(INB,INB)*(VENT(INB,INB)-V(INB))/ |
|---|
| 1073 | 1 (PH(INB)-PH(INB+1)) |
|---|
| 1074 | FV(INB)=FV(INB)-DX |
|---|
| 1075 | FV(INB-1)=FV(INB-1)+DX*(PH(INB)-PH(INB+1))/ |
|---|
| 1076 | 1 (PH(INB-1)-PH(INB)) |
|---|
| 1077 | DO J=1,NTRA |
|---|
| 1078 | EX=0.1*MENT(INB,INB)*(TRAENT(INB,INB,J) |
|---|
| 1079 | 1 -TRA(INB,J))/(PH(INB)-PH(INB+1)) |
|---|
| 1080 | FTRA(INB,J)=FTRA(INB,J)-EX |
|---|
| 1081 | FTRA(INB-1,J)=FTRA(INB-1,J)+EX* |
|---|
| 1082 | 1 (PH(INB)-PH(INB+1))/(PH(INB-1)-PH(INB)) |
|---|
| 1083 | ENDDO |
|---|
| 1084 | C |
|---|
| 1085 | C *** HOMOGINIZE TENDENCIES BELOW CLOUD BASE *** |
|---|
| 1086 | C |
|---|
| 1087 | ASUM=0.0 |
|---|
| 1088 | BSUM=0.0 |
|---|
| 1089 | CSUM=0.0 |
|---|
| 1090 | DSUM=0.0 |
|---|
| 1091 | DO 650 I=1,ICB-1 |
|---|
| 1092 | ASUM=ASUM+FT(I)*(PH(I)-PH(I+1)) |
|---|
| 1093 | BSUM=BSUM+FR(I)*(LV(I)+(CL-CPD)*(T(I)-T(1)))* |
|---|
| 1094 | 1 (PH(I)-PH(I+1)) |
|---|
| 1095 | CSUM=CSUM+(LV(I)+(CL-CPD)*(T(I)-T(1)))*(PH(I)-PH(I+1)) |
|---|
| 1096 | DSUM=DSUM+T(I)*(PH(I)-PH(I+1))/TH(I) |
|---|
| 1097 | 650 CONTINUE |
|---|
| 1098 | DO 700 I=1,ICB-1 |
|---|
| 1099 | FT(I)=ASUM*T(I)/(TH(I)*DSUM) |
|---|
| 1100 | FR(I)=BSUM/CSUM |
|---|
| 1101 | 700 CONTINUE |
|---|
| 1102 | C |
|---|
| 1103 | C *** RESET COUNTER AND RETURN *** |
|---|
| 1104 | C |
|---|
| 1105 | SIG(ND)=2.0 |
|---|
| 1106 | c |
|---|
| 1107 | c |
|---|
| 1108 | do i = 1, nd |
|---|
| 1109 | upwd(i) = 0.0 |
|---|
| 1110 | dnwd(i) = 0.0 |
|---|
| 1111 | c sb dnwd0(i) = - mp(i) |
|---|
| 1112 | enddo |
|---|
| 1113 | c |
|---|
| 1114 | do i = 1, nl |
|---|
| 1115 | dnwd0(i) = - mp(i) |
|---|
| 1116 | enddo |
|---|
| 1117 | do i = nl+1, nd |
|---|
| 1118 | dnwd0(i) = 0. |
|---|
| 1119 | enddo |
|---|
| 1120 | c |
|---|
| 1121 | do i = icb, inb |
|---|
| 1122 | upwd(i) = 0.0 |
|---|
| 1123 | dnwd(i) = 0.0 |
|---|
| 1124 | |
|---|
| 1125 | do k =i, inb |
|---|
| 1126 | up1=0.0 |
|---|
| 1127 | dn1=0.0 |
|---|
| 1128 | do n = 1, i-1 |
|---|
| 1129 | up1 = up1 + ment(n,k) |
|---|
| 1130 | dn1 = dn1 - ment(k,n) |
|---|
| 1131 | enddo |
|---|
| 1132 | upwd(i) = upwd(i)+ m(k) + up1 |
|---|
| 1133 | dnwd(i) = dnwd(i) + dn1 |
|---|
| 1134 | enddo |
|---|
| 1135 | enddo |
|---|
| 1136 | |
|---|
| 1137 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1138 | c DETERMINATION DE LA VARIATION DE FLUX ASCENDANT ENTRE |
|---|
| 1139 | C DEUX NIVEAU NON DILUE Mike |
|---|
| 1140 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1141 | |
|---|
| 1142 | |
|---|
| 1143 | c sb do i=1,ND |
|---|
| 1144 | c sb Mike(i)=M(i) |
|---|
| 1145 | c sb enddo |
|---|
| 1146 | |
|---|
| 1147 | do i = 1, NL |
|---|
| 1148 | Mike(i) = M(i) |
|---|
| 1149 | enddo |
|---|
| 1150 | do i = NL+1, ND |
|---|
| 1151 | Mike(i) = 0. |
|---|
| 1152 | enddo |
|---|
| 1153 | |
|---|
| 1154 | do i=1,nd |
|---|
| 1155 | Ma(i)=0 |
|---|
| 1156 | enddo |
|---|
| 1157 | |
|---|
| 1158 | c sb do i=1,nd |
|---|
| 1159 | c sb do j=i,nd |
|---|
| 1160 | c sb Ma(i)=Ma(i)+M(j) |
|---|
| 1161 | c sb enddo |
|---|
| 1162 | c sb enddo |
|---|
| 1163 | |
|---|
| 1164 | do i = 1, NL |
|---|
| 1165 | do j = i, NL |
|---|
| 1166 | Ma(i) = Ma(i) + M(j) |
|---|
| 1167 | enddo |
|---|
| 1168 | enddo |
|---|
| 1169 | c |
|---|
| 1170 | do i = NL+1, ND |
|---|
| 1171 | Ma(i) = 0. |
|---|
| 1172 | enddo |
|---|
| 1173 | c |
|---|
| 1174 | do i=1,ICB-1 |
|---|
| 1175 | Ma(i)=0 |
|---|
| 1176 | enddo |
|---|
| 1177 | |
|---|
| 1178 | |
|---|
| 1179 | |
|---|
| 1180 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1181 | C ICB REPRESENTE DE NIVEAU OU SE TROUVE LA |
|---|
| 1182 | c BASE DU NUAGE , ET INB LE TOP DU NUAGE |
|---|
| 1183 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1184 | |
|---|
| 1185 | |
|---|
| 1186 | do i=1,ND |
|---|
| 1187 | Mke(i)=upwd(i)+dnwd(i) |
|---|
| 1188 | enddo |
|---|
| 1189 | c$$$cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1190 | c$$$ call writeg1d(1,klev,ma,'ma ','ma ') |
|---|
| 1191 | c$$$ call writeg1d(1,klev,upwd,'upwd ','upwd ') |
|---|
| 1192 | c$$$ call writeg1d(1,klev,dnwd,'dnwd ','dnwd ') |
|---|
| 1193 | c$$$ call writeg1d(1,klev,dnwd0,'dnwd0 ','dnwd0 ') |
|---|
| 1194 | c$$$ call writeg1d(1,klev,tvp,'tvp ','tvp ') |
|---|
| 1195 | c$$$ call writeg1d(1,klev,tra(1:klev,3),'tra3 ','tra3 ') |
|---|
| 1196 | c$$$ call writeg1d(1,klev,tra(1:klev,4),'tra4 ','tra4 ') |
|---|
| 1197 | c$$$ call writeg1d(1,klev,tra(1:klev,5),'tra5 ','tra5 ') |
|---|
| 1198 | c$$$ call writeg1d(1,klev,tra(1:klev,6),'tra6 ','tra6 ') |
|---|
| 1199 | c$$$ call writeg1d(1,klev,tra(1:klev,7),'tra7 ','tra7 ') |
|---|
| 1200 | c$$$ call writeg1d(1,klev,tra(1:klev,8),'tra8 ','tra8 ') |
|---|
| 1201 | c$$$ call writeg1d(1,klev,tra(1:klev,9),'tra9 ','tra9 ') |
|---|
| 1202 | c$$$ call writeg1d(1,klev,tra(1:klev,10),'tra10','tra10') |
|---|
| 1203 | c$$$ call writeg1d(1,klev,tra(1:klev,11),'tra11','tra11') |
|---|
| 1204 | c$$$ call writeg1d(1,klev,tra(1:klev,12),'tra12','tra12') |
|---|
| 1205 | c$$$ call writeg1d(1,klev,tra(1:klev,13),'tra13','tra13') |
|---|
| 1206 | c$$$ call writeg1d(1,klev,tra(1:klev,14),'tra14','tra14') |
|---|
| 1207 | c$$$ call writeg1d(1,klev,tra(1:klev,15),'tra15','tra15') |
|---|
| 1208 | c$$$ call writeg1d(1,klev,tra(1:klev,16),'tra16','tra16') |
|---|
| 1209 | c$$$ call writeg1d(1,klev,tra(1:klev,17),'tra17','tra17') |
|---|
| 1210 | c$$$ call writeg1d(1,klev,tra(1:klev,18),'tra18','tra18') |
|---|
| 1211 | c$$$ call writeg1d(1,klev,tra(1:klev,19),'tra19','tra19') |
|---|
| 1212 | c$$$ call writeg1d(1,klev,tra(1:klev,20),'tra20','tra20 ') |
|---|
| 1213 | c$$$ call writeg1d(1,klev,trap(1:klev,1),'trp1','trp1') |
|---|
| 1214 | c$$$ call writeg1d(1,klev,trap(1:klev,2),'trp2','trp2') |
|---|
| 1215 | c$$$ call writeg1d(1,klev,trap(1:klev,3),'trp3','trp3') |
|---|
| 1216 | c$$$ call writeg1d(1,klev,trap(1:klev,4),'trp4','trp4') |
|---|
| 1217 | c$$$ call writeg1d(1,klev,trap(1:klev,5),'trp5','trp5') |
|---|
| 1218 | c$$$ call writeg1d(1,klev,trap(1:klev,10),'trp10','trp10') |
|---|
| 1219 | c$$$ call writeg1d(1,klev,trap(1:klev,12),'trp12','trp12') |
|---|
| 1220 | c$$$ call writeg1d(1,klev,trap(1:klev,15),'trp15','trp15') |
|---|
| 1221 | c$$$ call writeg1d(1,klev,trap(1:klev,20),'trp20','trp20') |
|---|
| 1222 | c$$$ call writeg1d(1,klev,ftra(1:klev,1),'ftr1 ','ftr1 ') |
|---|
| 1223 | c$$$ call writeg1d(1,klev,ftra(1:klev,2),'ftr2 ','ftr2 ') |
|---|
| 1224 | c$$$ call writeg1d(1,klev,ftra(1:klev,3),'ftr3 ','ftr3 ') |
|---|
| 1225 | c$$$ call writeg1d(1,klev,ftra(1:klev,4),'ftr4 ','ftr4 ') |
|---|
| 1226 | c$$$ call writeg1d(1,klev,ftra(1:klev,5),'ftr5 ','ftr5 ') |
|---|
| 1227 | c$$$ call writeg1d(1,klev,ftra(1:klev,6),'ftr6 ','ftr6 ') |
|---|
| 1228 | c$$$ call writeg1d(1,klev,ftra(1:klev,7),'ftr7 ','ftr7 ') |
|---|
| 1229 | c$$$ call writeg1d(1,klev,ftra(1:klev,8),'ftr8 ','ftr8 ') |
|---|
| 1230 | c$$$ call writeg1d(1,klev,ftra(1:klev,9),'ftr9 ','ftr9 ') |
|---|
| 1231 | c$$$ call writeg1d(1,klev,ftra(1:klev,10),'ftr10','ftr10') |
|---|
| 1232 | c$$$ call writeg1d(1,klev,ftra(1:klev,11),'ftr11','ftr11') |
|---|
| 1233 | c$$$ call writeg1d(1,klev,ftra(1:klev,12),'ftr12','ftr12') |
|---|
| 1234 | c$$$ call writeg1d(1,klev,ftra(1:klev,13),'ftr13','ftr13') |
|---|
| 1235 | c$$$ call writeg1d(1,klev,ftra(1:klev,14),'ftr14','ftr14') |
|---|
| 1236 | c$$$ call writeg1d(1,klev,ftra(1:klev,15),'ftr15','ftr15') |
|---|
| 1237 | c$$$ call writeg1d(1,klev,ftra(1:klev,16),'ftr16','ftr16') |
|---|
| 1238 | c$$$ call writeg1d(1,klev,ftra(1:klev,17),'ftr17','ftr17') |
|---|
| 1239 | c$$$ call writeg1d(1,klev,ftra(1:klev,18),'ftr18','ftr18') |
|---|
| 1240 | c$$$ call writeg1d(1,klev,ftra(1:klev,19),'ftr19','ftr19') |
|---|
| 1241 | c$$$ call writeg1d(1,klev,ftra(1:klev,20),'ftr20','ftr20 ') |
|---|
| 1242 | c$$$ call writeg1d(1,klev,mp,'mp ','mp ') |
|---|
| 1243 | c$$$ call writeg1d(1,klev,Mke,'Mke ','Mke ') |
|---|
| 1244 | |
|---|
| 1245 | |
|---|
| 1246 | |
|---|
| 1247 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1248 | c |
|---|
| 1249 | c |
|---|
| 1250 | RETURN |
|---|
| 1251 | END |
|---|
| 1252 | C --------------------------------------------------------------------------- |
|---|
