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