| 1 |
C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.1 2004/06/15 14:47:23 molod Exp $ |
| 2 |
C $Name: $ |
| 3 |
|
| 4 |
#include "CPP_OPTIONS.h" |
| 5 |
subroutine moistio (ndmoist,istrip,npcs,pz,tz,qz,ntracer,ptracer, |
| 6 |
. pkht,qqz,dumoist,dvmoist,dtmoist,dqmoist, |
| 7 |
. im,jm,lm,sige,sig,dsig,ptop, |
| 8 |
. iras,rainlsp,rainconv,snowfall, |
| 9 |
. nswcld,cldtot_sw,cldras_sw,cldlsp_sw,nswlz,swlz, |
| 10 |
. nlwcld,cldtot_lw,cldras_lw,cldlsp_lw,nlwlz,lwlz, |
| 11 |
. lpnt,myid) |
| 12 |
|
| 13 |
#ifdef ALLOW_DIAGNOSTICS |
| 14 |
#include "diagnostics.h" |
| 15 |
#endif |
| 16 |
|
| 17 |
c Input Variables |
| 18 |
c --------------- |
| 19 |
integer ndmoist,istrip,npcs,myid |
| 20 |
|
| 21 |
integer im,jm,lm |
| 22 |
real ptop |
| 23 |
real sige(lm+1) |
| 24 |
real sig(lm) |
| 25 |
real dsig(lm) |
| 26 |
|
| 27 |
integer ntracer,ptracer |
| 28 |
|
| 29 |
real pz(im,jm) |
| 30 |
real tz(im,jm,lm) |
| 31 |
real qz(im,jm,lm,ntracer) |
| 32 |
|
| 33 |
real pkht(im,jm,lm) |
| 34 |
|
| 35 |
real qqz(im,jm,lm) |
| 36 |
|
| 37 |
real dumoist(im,jm,lm) |
| 38 |
real dvmoist(im,jm,lm) |
| 39 |
real dtmoist(im,jm,lm) |
| 40 |
real dqmoist(im,jm,lm,ntracer) |
| 41 |
|
| 42 |
integer iras |
| 43 |
real rainlsp(im,jm) |
| 44 |
real rainconv(im,jm) |
| 45 |
real snowfall(im,jm) |
| 46 |
|
| 47 |
integer nswcld,nswlz |
| 48 |
real cldlsp_sw(im,jm,lm) |
| 49 |
real cldras_sw(im,jm,lm) |
| 50 |
real cldtot_sw(im,jm,lm) |
| 51 |
real swlz(im,jm,lm) |
| 52 |
|
| 53 |
integer nlwcld,nlwlz |
| 54 |
real cldlsp_lw(im,jm,lm) |
| 55 |
real cldras_lw(im,jm,lm) |
| 56 |
real cldtot_lw(im,jm,lm) |
| 57 |
real lwlz(im,jm,lm) |
| 58 |
|
| 59 |
logical lpnt |
| 60 |
|
| 61 |
c Local Variables |
| 62 |
c --------------- |
| 63 |
integer ncrnd,nsecf,nsubmax |
| 64 |
|
| 65 |
real fracqq, rh,temp1,temp2,dum |
| 66 |
integer snowcrit, lup |
| 67 |
parameter (fracqq = 0.1) |
| 68 |
|
| 69 |
real cldsr(im,jm,lm) |
| 70 |
real srcld(istrip,lm) |
| 71 |
|
| 72 |
real plev |
| 73 |
real cldnow,cldlsp_mem,cldras_mem,cldras,watnow,watmin,cldmin |
| 74 |
real cldprs(im,jm),cldtmp(im,jm) |
| 75 |
real cldhi (im,jm),cldlow(im,jm) |
| 76 |
real cldmid(im,jm),totcld(im,jm) |
| 77 |
integer midlevel,lowlevel |
| 78 |
|
| 79 |
real CLDLS(im,jm,lm) , CPEN(im,jm,lm) |
| 80 |
real tmpimjm(im,jm) |
| 81 |
real lsp_new(im,jm) |
| 82 |
real conv_new(im,jm) |
| 83 |
real snow_new(im,jm) |
| 84 |
|
| 85 |
real qqcolmin(im,jm) |
| 86 |
real qqcolmax(im,jm) |
| 87 |
integer levpbl(im,jm) |
| 88 |
|
| 89 |
c Gathered Arrays for Variable Cloud Base |
| 90 |
c --------------------------------------- |
| 91 |
real raincgath(im*jm) |
| 92 |
real pigather(im*jm) |
| 93 |
real thgather(im*jm,lm) |
| 94 |
real shgather(im*jm,lm) |
| 95 |
real pkzgather(im*jm,lm) |
| 96 |
real pkegather(im*jm,lm) |
| 97 |
real tmpgather(im*jm,lm) |
| 98 |
real deltgather(im*jm,lm) |
| 99 |
real delqgather(im*jm,lm) |
| 100 |
real ugather(im*jm,lm,ntracer) |
| 101 |
real delugather(im*jm,lm,ntracer) |
| 102 |
real deltrnev(im*jm,lm) |
| 103 |
real delqrnev(im*jm,lm) |
| 104 |
|
| 105 |
integer nindeces(lm) |
| 106 |
integer pblindex(im*jm) |
| 107 |
integer levgather(im*jm) |
| 108 |
|
| 109 |
c Stripped Arrays |
| 110 |
c --------------- |
| 111 |
real saveth (istrip,lm) |
| 112 |
real saveq (istrip,lm) |
| 113 |
real saveu (istrip,lm,ntracer) |
| 114 |
real usubcl (istrip, ntracer) |
| 115 |
|
| 116 |
real ple(istrip,lm+1), gam(istrip,lm) |
| 117 |
real TL(ISTRIP,lm) , SHL(ISTRIP,lm) |
| 118 |
real PL(ISTRIP,lm) , PLK(ISTRIP,lm) |
| 119 |
real PLKE(ISTRIP,lm+1) |
| 120 |
real TH(ISTRIP,lm) ,CVTH(ISTRIP,lm) |
| 121 |
real SHSAT(ISTRIP,lm) , CVQ(ISTRIP,lm) |
| 122 |
real UL(ISTRIP,lm,ntracer) |
| 123 |
real cvu(istrip,lm,ntracer) |
| 124 |
real CLMAXO(ISTRIP,lm),CLBOTH(ISTRIP,lm) |
| 125 |
real CLSBTH(ISTRIP,lm) |
| 126 |
real TMP1(ISTRIP,lm), TMP2(ISTRIP,lm) |
| 127 |
real TMP3(ISTRIP,lm), TMP4(ISTRIP,lm+1) |
| 128 |
real TMP5(ISTRIP,lm+1) |
| 129 |
integer ITMP1(ISTRIP,lm), ITMP2(ISTRIP,lm) |
| 130 |
integer ITMP3(ISTRIP,lm) |
| 131 |
|
| 132 |
real PRECIP(ISTRIP), PCMID(ISTRIP), PCNET(ISTRIP) |
| 133 |
real PCLOW (ISTRIP), SP(ISTRIP), PREP(ISTRIP) |
| 134 |
real PCPEN (ISTRIP,lm) |
| 135 |
integer pbl(istrip),depths(lm) |
| 136 |
|
| 137 |
real cldlz(istrip,lm), cldwater(im,jm,lm) |
| 138 |
real rhfrac(istrip), rhmin, pup, ppbl, rhcrit(istrip,lm) |
| 139 |
real offset, alpha, rasmax |
| 140 |
|
| 141 |
logical first |
| 142 |
logical lras |
| 143 |
real clfrac (istrip,lm) |
| 144 |
real cldmas (istrip,lm) |
| 145 |
real detrain(istrip,lm) |
| 146 |
real psubcld (istrip), psubcldg (im,jm) |
| 147 |
real psubcld_cnt(istrip), psubcldgc(im,jm) |
| 148 |
real rnd(lm/2) |
| 149 |
DATA FIRST /.TRUE./ |
| 150 |
|
| 151 |
integer imstp,nltop,nsubcl,nlras,nsubmin |
| 152 |
integer i,j,iloop,index,l,nn,num,numdeps,nt |
| 153 |
real tmstp,tminv,sday,grav,alhl,cp,elocp,gamfac |
| 154 |
real rkappa,p0kappa,p0kinv,ptopkap,pcheck |
| 155 |
real tice,getcon,pi |
| 156 |
|
| 157 |
C ********************************************************************** |
| 158 |
C **** INITIALIZATION **** |
| 159 |
C ********************************************************************** |
| 160 |
|
| 161 |
IMSTP = nsecf(NDMOIST) |
| 162 |
TMSTP = FLOAT(IMSTP) |
| 163 |
TMINV = 1. / TMSTP |
| 164 |
|
| 165 |
C Minimum Large-Scale Cloud Fraction at rhcrit |
| 166 |
alpha = 0.80 |
| 167 |
C Difference in fraction between SR and LS Threshold |
| 168 |
offset = 0.10 |
| 169 |
C Large-Scale Relative Humidity Threshold in PBL |
| 170 |
rhmin = 0.90 |
| 171 |
C Maximum Cloud Fraction associated with RAS |
| 172 |
rasmax = 1.00 |
| 173 |
|
| 174 |
nn = 3*3600.0/tmstp + 1 |
| 175 |
C Threshold for Cloud Fraction Memory |
| 176 |
cldmin = rasmax*(1.0-tmstp/3600.)**nn |
| 177 |
C Threshold for Cloud Liquid Water Memory |
| 178 |
watmin = 1.0e-8 |
| 179 |
|
| 180 |
SDAY = GETCON('SDAY') |
| 181 |
GRAV = GETCON('GRAVITY') |
| 182 |
ALHL = GETCON('LATENT HEAT COND') |
| 183 |
CP = GETCON('CP') |
| 184 |
ELOCP = GETCON('LATENT HEAT COND') / GETCON('CP') |
| 185 |
GAMFAC = GETCON('LATENT HEAT COND') * GETCON('EPS') * ELOCP |
| 186 |
. / GETCON('RGAS') |
| 187 |
RKAPPA = GETCON('KAPPA') |
| 188 |
P0KAPPA = 1000.0**RKAPPA |
| 189 |
P0KINV = 1. / P0KAPPA |
| 190 |
PTOPKAP = PTOP**RKAPPA |
| 191 |
tice = getcon('FREEZING-POINT') |
| 192 |
PI = 4.*atan(1.) |
| 193 |
|
| 194 |
c Determine Upper Level for Cumulus Convection |
| 195 |
c and Total number of Random Clouds to Check |
| 196 |
c --------------------------------------------- |
| 197 |
|
| 198 |
NLTOP = 1 |
| 199 |
DO L=1,lm |
| 200 |
PCHECK = (1000.-ptop)*SIG(L) + PTOP |
| 201 |
IF (PCHECK.GE.10.) THEN |
| 202 |
NLTOP = L |
| 203 |
GO TO 2 |
| 204 |
ENDIF |
| 205 |
ENDDO |
| 206 |
2 CONTINUE |
| 207 |
ncrnd = (lm-nltop+1)/2 |
| 208 |
|
| 209 |
c Determine Minimum Number of Levels in Sub-Cloud (50 mb) Layer |
| 210 |
c ------------------------------------------------------------- |
| 211 |
nsubmin = lm |
| 212 |
nsubmax = 1 |
| 213 |
DO L=lm-1,1,-1 |
| 214 |
PCHECK = (1000.-ptop)*SIG(L) + PTOP |
| 215 |
IF( PCHECK.GE.950.0 ) nsubmin = L |
| 216 |
IF( PCHECK.GE.750.0 ) nsubmax = L |
| 217 |
ENDDO |
| 218 |
|
| 219 |
if(first .and. myid.eq.0) then |
| 220 |
print * |
| 221 |
print *,'Top Level Allowed for Convection : ',nltop, |
| 222 |
. ' (',(1000.-ptop)*SIG(nltop) + PTOP,' mb)' |
| 223 |
print *,' Highest Sub-Cloud Level: ',nsubmax, |
| 224 |
. ' (',(1000.-ptop)*SIG(nsubmax) + PTOP,' mb)' |
| 225 |
print *,' Lowest Sub-Cloud Level: ',nsubmin, |
| 226 |
. ' (',(1000.-ptop)*SIG(nsubmin) + PTOP,' mb)' |
| 227 |
print *,' Total Number of Random Clouds: ',ncrnd |
| 228 |
print * |
| 229 |
first = .false. |
| 230 |
endif |
| 231 |
|
| 232 |
c And now find PBL depth - the level where qq = fracqq * qq at surface |
| 233 |
c -------------------------------------------------------------------- |
| 234 |
do j = 1,jm |
| 235 |
do i = 1,im |
| 236 |
qqcolmin(i,j) = qqz(i,j,lm)*fracqq |
| 237 |
qqcolmax(i,j) = qqz(i,j,lm) |
| 238 |
levpbl(i,j) = lm+1 |
| 239 |
enddo |
| 240 |
enddo |
| 241 |
|
| 242 |
DO L = lm-1,1,-1 |
| 243 |
DO j = 1,jm |
| 244 |
DO i = 1,im |
| 245 |
IF((qqz(i,j,l).gt.qqcolmax(i,j)) |
| 246 |
1 .and.(levpbl(i,j).eq.lm+1))then |
| 247 |
qqcolmax(i,j) = qqz(i,j,l) |
| 248 |
qqcolmin(i,j) = fracqq*qqcolmax(i,j) |
| 249 |
endif |
| 250 |
if((qqz(i,j,l).lt.qqcolmin(i,j)) |
| 251 |
1 .and.(levpbl(i,j).eq.lm+1))levpbl(i,j)=L+1 |
| 252 |
enddo |
| 253 |
enddo |
| 254 |
enddo |
| 255 |
|
| 256 |
do j = 1,jm |
| 257 |
do i = 1,im |
| 258 |
if(levpbl(i,j).gt.nsubmin) levpbl(i,j) = nsubmin |
| 259 |
if(levpbl(i,j).lt.nsubmax) levpbl(i,j) = nsubmax |
| 260 |
enddo |
| 261 |
enddo |
| 262 |
|
| 263 |
|
| 264 |
c Set up the array of indeces of subcloud levels for the gathering |
| 265 |
c ---------------------------------------------------------------- |
| 266 |
index = 0 |
| 267 |
do L = nsubmin,nltop,-1 |
| 268 |
do j = 1,jm |
| 269 |
do i = 1,im |
| 270 |
if(levpbl(i,j).eq.L) then |
| 271 |
index = index + 1 |
| 272 |
pblindex(index) = (j-1)*im + i |
| 273 |
endif |
| 274 |
enddo |
| 275 |
enddo |
| 276 |
enddo |
| 277 |
|
| 278 |
do index = 1,im*jm |
| 279 |
levgather(index) = levpbl(pblindex(index),1) |
| 280 |
pigather(index) = pz(pblindex(index),1) |
| 281 |
enddo |
| 282 |
|
| 283 |
do L = 1,lm |
| 284 |
do index = 1,im*jm |
| 285 |
thgather(index,L) = tz(pblindex(index),1,L) |
| 286 |
shgather(index,L) = qz(pblindex(index),1,L,1) |
| 287 |
pkegather(index,L) = pkht(pblindex(index),1,L) |
| 288 |
enddo |
| 289 |
enddo |
| 290 |
do nt = 1,ntracer-ptracer |
| 291 |
do L = 1,lm |
| 292 |
do index = 1,im*jm |
| 293 |
ugather(index,L,nt) = qz(pblindex(index),1,L,nt+ptracer) |
| 294 |
enddo |
| 295 |
enddo |
| 296 |
enddo |
| 297 |
|
| 298 |
call pkappa(pigather,pkegather,pkzgather,ptop,sige,dsig,im,jm,lm) |
| 299 |
|
| 300 |
c bump the counter for number of calls to convection |
| 301 |
c -------------------------------------------------- |
| 302 |
iras = iras + 1 |
| 303 |
if( iras.ge.1e9 ) iras = 1 |
| 304 |
|
| 305 |
c select the 'random' cloud detrainment levels for RAS |
| 306 |
c ---------------------------------------------------- |
| 307 |
call rndcloud(iras,ncrnd,rnd,myid) |
| 308 |
|
| 309 |
do l=1,lm |
| 310 |
do j=1,jm |
| 311 |
do i=1,im |
| 312 |
dtmoist(i,j,l) = 0. |
| 313 |
do nt = 1,ntracer |
| 314 |
dqmoist(i,j,l,nt) = 0. |
| 315 |
enddo |
| 316 |
enddo |
| 317 |
enddo |
| 318 |
enddo |
| 319 |
|
| 320 |
C*********************************************************************** |
| 321 |
C **** LOOP OVER NPCS PEICES **** |
| 322 |
C ********************************************************************** |
| 323 |
|
| 324 |
DO 1000 NN = 1,NPCS |
| 325 |
|
| 326 |
C ********************************************************************** |
| 327 |
C **** VARIABLE INITIALIZATION **** |
| 328 |
C ********************************************************************** |
| 329 |
|
| 330 |
CALL STRIP ( pigather, SP ,im*jm,ISTRIP,1 ,NN ) |
| 331 |
CALL STRIP ( pkzgather, PLK ,im*jm,ISTRIP,lm,NN ) |
| 332 |
CALL STRIP ( pkegather, PLKE ,im*jm,ISTRIP,lm,NN ) |
| 333 |
CALL STRIP ( thgather, TH ,im*jm,ISTRIP,lm,NN ) |
| 334 |
CALL STRIP ( shgather, SHL ,im*jm,ISTRIP,lm,NN ) |
| 335 |
CALL STRINT( levgather, pbl ,im*jm,ISTRIP,1 ,NN ) |
| 336 |
|
| 337 |
do nt = 1,ntracer-ptracer |
| 338 |
call strip ( ugather(1,1,nt), ul(1,1,nt),im*jm,istrip,lm,nn ) |
| 339 |
enddo |
| 340 |
|
| 341 |
do l = 1,lm |
| 342 |
do i = 1,istrip |
| 343 |
PL(I,L) = SIG(L)*SP(I) + PTOP |
| 344 |
PLE(I,L) = SIGE(L)*SP(I) + PTOP |
| 345 |
enddo |
| 346 |
enddo |
| 347 |
|
| 348 |
do i = 1,istrip |
| 349 |
PLE(I,lm+1) = SP(I) + PTOP |
| 350 |
enddo |
| 351 |
|
| 352 |
C ********************************************************************** |
| 353 |
C **** SETUP FOR RAS CUMULUS PARAMETERIZATION **** |
| 354 |
C ********************************************************************** |
| 355 |
|
| 356 |
DO L = 1,lm |
| 357 |
DO I = 1,ISTRIP |
| 358 |
TH(I,L) = TH(I,L) * P0KAPPA |
| 359 |
CLMAXO(I,L) = 0. |
| 360 |
CLBOTH(I,L) = 0. |
| 361 |
cldmas(I,L) = 0. |
| 362 |
detrain(I,L) = 0. |
| 363 |
ENDDO |
| 364 |
ENDDO |
| 365 |
|
| 366 |
do L = 1,lm |
| 367 |
depths(L) = 0 |
| 368 |
enddo |
| 369 |
|
| 370 |
numdeps = 0 |
| 371 |
do L = nsubmin,nltop,-1 |
| 372 |
nindeces(L) = 0 |
| 373 |
do i = 1,istrip |
| 374 |
if(pbl(i).eq.L) nindeces(L) = nindeces(L) + 1 |
| 375 |
enddo |
| 376 |
if(nindeces(L).gt.0) then |
| 377 |
numdeps = numdeps + 1 |
| 378 |
depths(numdeps) = L |
| 379 |
endif |
| 380 |
enddo |
| 381 |
|
| 382 |
|
| 383 |
C Initiate a do-loop around RAS for the number of different |
| 384 |
C sub-cloud layer depths in this strip |
| 385 |
C --If all subcloud depths are the same, execute loop once |
| 386 |
C Otherwise loop over different subcloud layer depths |
| 387 |
|
| 388 |
num = 1 |
| 389 |
DO iloop = 1,numdeps |
| 390 |
|
| 391 |
nsubcl = depths(iloop) |
| 392 |
|
| 393 |
c Compute sub-cloud values for Temperature and Spec.Hum. |
| 394 |
c ------------------------------------------------------ |
| 395 |
DO 600 I=num,num+nindeces(nsubcl)-1 |
| 396 |
TMP1(I,2) = 0. |
| 397 |
TMP1(I,3) = 0. |
| 398 |
600 CONTINUE |
| 399 |
|
| 400 |
NLRAS = NSUBCL - NLTOP + 1 |
| 401 |
DO 601 L=NSUBCL,lm |
| 402 |
DO 602 I=num,num+nindeces(nsubcl)-1 |
| 403 |
TMP1(I,2) = TMP1(I,2) + (PLE(I,L+1)-PLE(I,L))*TH (I,L)/sp(i) |
| 404 |
TMP1(I,3) = TMP1(I,3) + (PLE(I,L+1)-PLE(I,L))*SHL(I,L)/sp(i) |
| 405 |
602 CONTINUE |
| 406 |
601 CONTINUE |
| 407 |
DO 603 I=num,num+nindeces(nsubcl)-1 |
| 408 |
TMP1(I,4) = 1. / ( (PLE(I,lm+1)-PLE(I,NSUBCL))/sp(I) ) |
| 409 |
TH(I,NSUBCL) = TMP1(I,2)*TMP1(I,4) |
| 410 |
SHL(I,NSUBCL) = TMP1(I,3)*TMP1(I,4) |
| 411 |
603 CONTINUE |
| 412 |
|
| 413 |
c Save initial value of tracers and compute sub-cloud value |
| 414 |
c --------------------------------------------------------- |
| 415 |
DO NT = 1,ntracer-ptracer |
| 416 |
do L = 1,lm |
| 417 |
do i = num,num+nindeces(nsubcl)-1 |
| 418 |
saveu(i,L,nt) = ul(i,L,nt) |
| 419 |
enddo |
| 420 |
enddo |
| 421 |
DO I=num,num+nindeces(nsubcl)-1 |
| 422 |
TMP1(I,2) = 0. |
| 423 |
ENDDO |
| 424 |
DO L=NSUBCL,lm |
| 425 |
DO I=num,num+nindeces(nsubcl)-1 |
| 426 |
TMP1(I,2) = TMP1(I,2)+(PLE(I,L+1)-PLE(I,L))*UL(I,L,NT)/sp(i) |
| 427 |
ENDDO |
| 428 |
ENDDO |
| 429 |
DO I=num,num+nindeces(nsubcl)-1 |
| 430 |
UL(I,NSUBCL,NT) = TMP1(I,2)*TMP1(I,4) |
| 431 |
usubcl(i,nt) = ul(i,nsubcl,nt) |
| 432 |
ENDDO |
| 433 |
ENDDO |
| 434 |
|
| 435 |
c Compute Pressure Arrays for RAS |
| 436 |
c ------------------------------- |
| 437 |
DO 111 L=1,lm |
| 438 |
DO 112 I=num,num+nindeces(nsubcl)-1 |
| 439 |
TMP4(I,L) = PLE(I,L) |
| 440 |
112 CONTINUE |
| 441 |
111 CONTINUE |
| 442 |
DO I=num,num+nindeces(nsubcl)-1 |
| 443 |
TMP5(I,1) = PTOPKAP / P0KAPPA |
| 444 |
ENDDO |
| 445 |
DO L=2,lm |
| 446 |
DO I=num,num+nindeces(nsubcl)-1 |
| 447 |
TMP5(I,L) = PLKE(I,L-1)*P0KINV |
| 448 |
ENDDO |
| 449 |
ENDDO |
| 450 |
DO I=num,num+nindeces(nsubcl)-1 |
| 451 |
TMP4(I,lm+1) = PLE (I,lm+1) |
| 452 |
TMP5(I,lm+1) = PLKE(I,lm)*P0KINV |
| 453 |
ENDDO |
| 454 |
DO 113 I=num,num+nindeces(nsubcl)-1 |
| 455 |
TMP4(I,NSUBCL+1) = PLE (I,lm+1) |
| 456 |
TMP5(I,NSUBCL+1) = PLKE(I,lm)*P0KINV |
| 457 |
113 CONTINUE |
| 458 |
|
| 459 |
do i=num,num+nindeces(nsubcl)-1 |
| 460 |
C Temperature at top of sub-cloud layer |
| 461 |
tmp2(i,1) = TH(i,NSUBCL) * PLKE(i,NSUBCL)/P0KAPPA |
| 462 |
C Pressure at top of sub-cloud layer |
| 463 |
tmp2(i,2) = tmp4(i,nsubcl) |
| 464 |
enddo |
| 465 |
|
| 466 |
C CHANGED THIS: no RH requirement for RAS |
| 467 |
c call vqsat ( tmp2(num,1),tmp2(num,2),tmp2(num,3), |
| 468 |
c . dum,.false.,nindeces(nsubcl) ) |
| 469 |
c do i=num,num+nindeces(nsubcl)-1 |
| 470 |
c rh = SHL(I,NSUBCL) / tmp2(i,3) |
| 471 |
c if (rh .le. 0.85) then |
| 472 |
c rhfrac(i) = 0. |
| 473 |
c else if (rh .ge. 0.95) then |
| 474 |
c rhfrac(i) = 1. |
| 475 |
c else |
| 476 |
c rhfrac(i) = (rh-0.85)*10. |
| 477 |
c endif |
| 478 |
c enddo |
| 479 |
do i=num,num+nindeces(nsubcl)-1 |
| 480 |
rhfrac(i) = 1. |
| 481 |
enddo |
| 482 |
|
| 483 |
C Compute RH threshold for Large-scale condensation |
| 484 |
C Used in Slingo-Ritter clouds as well - define offset between SR and LS |
| 485 |
|
| 486 |
C Top level of atan func above this rh_threshold = rhmin |
| 487 |
pup = 600. |
| 488 |
do i=num,num+nindeces(nsubcl)-1 |
| 489 |
do L = nsubcl, lm |
| 490 |
rhcrit(i,L) = 1. |
| 491 |
enddo |
| 492 |
do L = 1, nsubcl-1 |
| 493 |
pcheck = (1000.-ptop)*sig(L) + ptop |
| 494 |
if (pcheck .le. pup) then |
| 495 |
rhcrit(i,L) = rhmin |
| 496 |
else |
| 497 |
ppbl = (1000.-ptop)*sig(nsubcl) + ptop |
| 498 |
rhcrit(i,L) = rhmin + (1.-rhmin)/(19.) * |
| 499 |
. ((atan( (2.*(pcheck-pup)/(ppbl-pup)-1.) * |
| 500 |
. tan(20.*pi/21.-0.5*pi) ) |
| 501 |
. + 0.5*pi) * 21./pi - 1.) |
| 502 |
endif |
| 503 |
enddo |
| 504 |
enddo |
| 505 |
|
| 506 |
c Save Initial Values of Temperature and Specific Humidity |
| 507 |
c -------------------------------------------------------- |
| 508 |
do L = 1,lm |
| 509 |
do i = num,num+nindeces(nsubcl)-1 |
| 510 |
saveth(i,L) = th (i,L) |
| 511 |
saveq (i,L) = shl(i,L) |
| 512 |
PCPEN (i,L) = 0. |
| 513 |
CLFRAC(i,L) = 0. |
| 514 |
enddo |
| 515 |
enddo |
| 516 |
|
| 517 |
CALL RAS ( NN,istrip,nindeces(nsubcl),NLRAS,NLTOP,lm,TMSTP |
| 518 |
1, UL(num,1,1),ntracer-ptracer,TH(num,NLTOP),SHL(num,NLTOP) |
| 519 |
2, TMP4(num,NLTOP), TMP5(num,NLTOP),rnd, ncrnd, PCPEN(num,NLTOP) |
| 520 |
3, CLBOTH(num,NLTOP), CLFRAC(num,NLTOP) |
| 521 |
4, cldmas(num,nltop), detrain(num,nltop) |
| 522 |
8, cp,grav,rkappa,alhl,rhfrac(num),rasmax ) |
| 523 |
|
| 524 |
c Compute Diagnostic CLDMAS in RAS Subcloud Layers |
| 525 |
c ------------------------------------------------ |
| 526 |
do L=nsubcl,lm |
| 527 |
dum = dsig(L)/(1.0-sige(nsubcl)) |
| 528 |
do I=num,num+nindeces(nsubcl)-1 |
| 529 |
cldmas(i,L) = cldmas(i,L-1) - dum*cldmas(i,nsubcl-1) |
| 530 |
enddo |
| 531 |
enddo |
| 532 |
|
| 533 |
c Update Theta and Moisture due to RAS |
| 534 |
c ------------------------------------ |
| 535 |
DO L=1,nsubcl |
| 536 |
DO I=num,num+nindeces(nsubcl)-1 |
| 537 |
CVTH(I,L) = (TH (I,L) - saveth(i,l)) |
| 538 |
CVQ (I,L) = (SHL(I,L) - saveq (i,l)) |
| 539 |
ENDDO |
| 540 |
ENDDO |
| 541 |
DO L=nsubcl+1,lm |
| 542 |
DO I=num,num+nindeces(nsubcl)-1 |
| 543 |
CVTH(I,L) = cvth(i,nsubcl) |
| 544 |
CVQ (I,L) = cvq (i,nsubcl) |
| 545 |
ENDDO |
| 546 |
ENDDO |
| 547 |
|
| 548 |
DO L=nsubcl+1,lm |
| 549 |
DO I=num,num+nindeces(nsubcl)-1 |
| 550 |
TH (I,L) = saveth(i,l) + cvth(i,l) |
| 551 |
SHL(I,L) = saveq (i,l) + cvq (i,l) |
| 552 |
ENDDO |
| 553 |
ENDDO |
| 554 |
DO L=1,lm |
| 555 |
DO I=num,num+nindeces(nsubcl)-1 |
| 556 |
CVTH(I,L) = CVTH(I,L) *P0KINV*SP(I)*tminv |
| 557 |
CVQ (I,L) = CVQ (I,L) *SP(I)*tminv |
| 558 |
ENDDO |
| 559 |
ENDDO |
| 560 |
|
| 561 |
c Compute Tracer Tendency due to RAS |
| 562 |
c ---------------------------------- |
| 563 |
do nt = 1,ntracer-ptracer |
| 564 |
DO L=1,nsubcl-1 |
| 565 |
DO I=num,num+nindeces(nsubcl)-1 |
| 566 |
CVU(I,L,nt) = ( UL(I,L,nt)-saveu(i,l,nt) )*sp(i)*tminv |
| 567 |
ENDDO |
| 568 |
ENDDO |
| 569 |
DO L=nsubcl,lm |
| 570 |
DO I=num,num+nindeces(nsubcl)-1 |
| 571 |
if( usubcl(i,nt).ne.0.0 ) then |
| 572 |
cvu(i,L,nt) = ( ul(i,nsubcl,nt)-usubcl(i,nt) ) * |
| 573 |
. ( saveu(i,L,nt)/usubcl(i,nt) )*sp(i)*tminv |
| 574 |
else |
| 575 |
cvu(i,L,nt) = 0.0 |
| 576 |
endif |
| 577 |
ENDDO |
| 578 |
ENDDO |
| 579 |
enddo |
| 580 |
|
| 581 |
c Compute Diagnostic PSUBCLD (Subcloud Layer Pressure) |
| 582 |
c ---------------------------------------------------- |
| 583 |
do i=num,num+nindeces(nsubcl)-1 |
| 584 |
lras = .false. |
| 585 |
do L=nltop,nsubcl |
| 586 |
if( cvq(i,L).ne.0.0 ) lras = .true. |
| 587 |
enddo |
| 588 |
psubcld (i) = 0.0 |
| 589 |
psubcld_cnt(i) = 0.0 |
| 590 |
if( lras ) then |
| 591 |
psubcld (i) = sp(i)+ptop-ple(i,nsubcl) |
| 592 |
psubcld_cnt(i) = 1.0 |
| 593 |
endif |
| 594 |
enddo |
| 595 |
|
| 596 |
|
| 597 |
C End of subcloud layer depth loop (iloop) |
| 598 |
|
| 599 |
num = num+nindeces(nsubcl) |
| 600 |
|
| 601 |
ENDDO |
| 602 |
|
| 603 |
C ********************************************************************** |
| 604 |
C **** TENDENCY UPDATES **** |
| 605 |
C **** (Keep 'Gathered' tendencies in 'gather' arrays now) **** |
| 606 |
C ********************************************************************** |
| 607 |
|
| 608 |
call paste( CVTH,deltgather,istrip,im*jm,lm,NN ) |
| 609 |
call paste( CVQ,delqgather,istrip,im*jm,lm,NN ) |
| 610 |
do nt = 1,ntracer-ptracer |
| 611 |
call paste( CVU(1,1,nt),delugather(1,1,nt),istrip,im*jm,lm,NN ) |
| 612 |
enddo |
| 613 |
|
| 614 |
C ********************************************************************** |
| 615 |
C And now paste some arrays for filling diagnostics |
| 616 |
C (use pkegather to hold detrainment and tmpgather for cloud mass flux) |
| 617 |
C ********************************************************************** |
| 618 |
|
| 619 |
if(icldmas .gt.0) call paste( cldmas,tmpgather,istrip,im*jm,lm,NN) |
| 620 |
if(idtrain .gt.0) call paste(detrain,pkegather,istrip,im*jm,lm,NN) |
| 621 |
if(ipsubcld.gt.0) then |
| 622 |
call paste(psubcld ,psubcldg ,istrip,im*jm,1,NN) |
| 623 |
call paste(psubcld_cnt,psubcldgc,istrip,im*jm,1,NN) |
| 624 |
endif |
| 625 |
|
| 626 |
C ********************************************************************* |
| 627 |
C **** RE-EVAPORATION OF PENETRATING CONVECTIVE RAIN **** |
| 628 |
C ********************************************************************* |
| 629 |
|
| 630 |
CALL STRIP ( thgather,TH ,im*jm,ISTRIP,lm,NN) |
| 631 |
CALL STRIP ( shgather,SHL,im*jm,ISTRIP,lm,NN) |
| 632 |
DO L=1,lm |
| 633 |
DO I=1,ISTRIP |
| 634 |
TH(I,L) = TH(I,L) + CVTH(I,L)*tmstp/SP(I) |
| 635 |
SHL(I,L) = SHL(I,L) + CVQ(I,L)*tmstp/SP(I) |
| 636 |
TL(I,L) = TH(I,L)*PLK(I,L) |
| 637 |
saveth(I,L) = th(I,L) |
| 638 |
saveq (I,L) = SHL(I,L) |
| 639 |
ENDDO |
| 640 |
ENDDO |
| 641 |
|
| 642 |
CALL RNEVP (NN,ISTRIP,lm,TL,SHL,PCPEN,PL,CLFRAC,SP,DSIG,PLKE, |
| 643 |
. PLK,TH,TMP1,TMP2,TMP3,ITMP1,ITMP2,PCNET,PRECIP, |
| 644 |
. CLSBTH,TMSTP,1.,cp,grav,alhl,gamfac,cldlz,rhcrit,offset,alpha) |
| 645 |
|
| 646 |
C ********************************************************************** |
| 647 |
C **** TENDENCY UPDATES **** |
| 648 |
C ********************************************************************** |
| 649 |
|
| 650 |
DO L=1,lm |
| 651 |
|
| 652 |
DO I =1,ISTRIP |
| 653 |
TMP1(I,L) = sp(i) * (SHL(I,L)-saveq(I,L)) * tminv |
| 654 |
ENDDO |
| 655 |
CALL PSTBMP(TMP1(1,L),delqgather(1,L),ISTRIP,im*jm,1,NN) |
| 656 |
|
| 657 |
DO I =1,ISTRIP |
| 658 |
TMP1(I,L) = sp(i) * ((TL(I,L)/PLK(I,L))-saveth(i,l)) * tminv |
| 659 |
ENDDO |
| 660 |
CALL PSTBMP(TMP1(1,L),deltgather(1,L),ISTRIP,im*jm,1,NN) |
| 661 |
|
| 662 |
C Paste rain evap tendencies into arrays for diagnostic output |
| 663 |
c ------------------------------------------------------------ |
| 664 |
if(idtls.gt.0)then |
| 665 |
DO I =1,ISTRIP |
| 666 |
TMP1(I,L) = ((TL(I,L)/PLK(I,L))-saveth(i,l))*plk(i,l)*sday*tminv |
| 667 |
ENDDO |
| 668 |
call paste(tmp1(1,L),deltrnev(1,L),istrip,im*jm,1,NN) |
| 669 |
endif |
| 670 |
|
| 671 |
if(idqls.gt.0)then |
| 672 |
DO I =1,ISTRIP |
| 673 |
TMP1(I,L) = (SHL(I,L)-saveq(I,L)) * 1000. * sday * tminv |
| 674 |
ENDDO |
| 675 |
call paste(tmp1(1,L),delqrnev(1,L),istrip,im*jm,1,NN) |
| 676 |
endif |
| 677 |
|
| 678 |
ENDDO |
| 679 |
|
| 680 |
C ********************************************************************* |
| 681 |
C Add Non-Precipitating Clouds where the relative |
| 682 |
C humidity is less than 100% |
| 683 |
C Apply Cloud Top Entrainment Instability |
| 684 |
C ********************************************************************* |
| 685 |
|
| 686 |
do L=1,lm |
| 687 |
do i=1,istrip |
| 688 |
srcld(i,L) = -clsbth(i,L) |
| 689 |
enddo |
| 690 |
enddo |
| 691 |
|
| 692 |
call srclouds (saveth,saveq,plk,pl,plke,clsbth,cldlz,istrip,lm, |
| 693 |
. rhcrit,offset,alpha) |
| 694 |
|
| 695 |
do L=1,lm |
| 696 |
do i=1,istrip |
| 697 |
srcld(i,L) = srcld(i,L)+clsbth(i,L) |
| 698 |
enddo |
| 699 |
enddo |
| 700 |
|
| 701 |
C ********************************************************************* |
| 702 |
C **** PASTE CLOUD AMOUNTS **** |
| 703 |
C ********************************************************************* |
| 704 |
|
| 705 |
call paste ( srcld, cldsr,istrip,im*jm,lm,nn ) |
| 706 |
call paste ( cldlz,cldwater,istrip,im*jm,lm,nn ) |
| 707 |
call paste ( clsbth, cldls,istrip,im*jm,lm,nn ) |
| 708 |
call paste ( clboth, cpen ,istrip,im*jm,lm,nn ) |
| 709 |
|
| 710 |
c compute Total Accumulated Precip for Landsurface Model |
| 711 |
c ------------------------------------------------------ |
| 712 |
do i = 1,istrip |
| 713 |
C Initialize Rainlsp, Rainconv and Snowfall |
| 714 |
tmp1(i,1) = 0.0 |
| 715 |
tmp1(i,2) = 0.0 |
| 716 |
tmp1(i,3) = 0.0 |
| 717 |
enddo |
| 718 |
|
| 719 |
do i = 1,istrip |
| 720 |
prep(i) = PRECIP(I) + PCNET(I) |
| 721 |
tmp1(i,1) = PRECIP(I) |
| 722 |
tmp1(i,2) = pcnet(i) |
| 723 |
enddo |
| 724 |
c |
| 725 |
c check whether there is snow |
| 726 |
c------------------------------------------------------- |
| 727 |
c snow algorthm: |
| 728 |
c if temperature profile from the surface level to 700 mb |
| 729 |
c uniformaly c below zero, then precipitation (total) is |
| 730 |
c snowfall. Else there is no snow. |
| 731 |
c For version of level 70, the sigma level corresponding |
| 732 |
c to 700mb (assume the surface pressure is 1000mb) is |
| 733 |
c the 13th level from the surface |
| 734 |
c Runhua Yang Aug. 24 98 |
| 735 |
c added pcheck for 700mb - sharon sept 18, 1998 |
| 736 |
c------------------------------------------------------- |
| 737 |
|
| 738 |
pup = 700. |
| 739 |
do L = lm, 1, -1 |
| 740 |
pcheck = (1000.-ptop)*sig(L) + ptop |
| 741 |
if (pcheck .ge. pup) then |
| 742 |
lup = L |
| 743 |
endif |
| 744 |
enddo |
| 745 |
do i = 1,istrip |
| 746 |
snowcrit=0 |
| 747 |
do l=lup,lm |
| 748 |
if (saveth(i,l)*plk(i,l).le. tice ) then |
| 749 |
snowcrit=snowcrit+1 |
| 750 |
endif |
| 751 |
enddo |
| 752 |
if (snowcrit .eq. (lm-lup+1)) then |
| 753 |
tmp1(i,3) = prep(i) |
| 754 |
tmp1(i,1)=0.0 |
| 755 |
tmp1(i,2)=0.0 |
| 756 |
endif |
| 757 |
enddo |
| 758 |
|
| 759 |
CALL paste (tmp1(1,1), lsp_new,ISTRIP,im*jm,1,NN) |
| 760 |
CALL paste (tmp1(1,2),conv_new,ISTRIP,im*jm,1,NN) |
| 761 |
CALL paste (tmp1(1,3),snow_new,ISTRIP,im*jm,1,NN) |
| 762 |
|
| 763 |
if(iprecon.gt.0) then |
| 764 |
CALL paste (pcnet,raincgath,ISTRIP,im*jm,1,NN) |
| 765 |
endif |
| 766 |
|
| 767 |
C ********************************************************************* |
| 768 |
C **** End Major Stripped Region **** |
| 769 |
C ********************************************************************* |
| 770 |
|
| 771 |
1000 CONTINUE |
| 772 |
|
| 773 |
C Large Scale Rainfall, Conv rain, and snowfall |
| 774 |
c --------------------------------------------- |
| 775 |
call back2grd ( lsp_new,pblindex, lsp_new,im*jm) |
| 776 |
call back2grd (conv_new,pblindex,conv_new,im*jm) |
| 777 |
call back2grd (snow_new,pblindex,snow_new,im*jm) |
| 778 |
|
| 779 |
if(iprecon.gt.0) then |
| 780 |
call back2grd (raincgath,pblindex,raincgath,im*jm) |
| 781 |
endif |
| 782 |
|
| 783 |
c Subcloud Layer Pressure |
| 784 |
c ----------------------- |
| 785 |
if(ipsubcld.gt.0) then |
| 786 |
call back2grd (psubcldg ,pblindex,psubcldg ,im*jm) |
| 787 |
call back2grd (psubcldgc,pblindex,psubcldgc,im*jm) |
| 788 |
endif |
| 789 |
|
| 790 |
do L = 1,lm |
| 791 |
C Delta theta,q, convective, max and ls clouds |
| 792 |
c -------------------------------------------- |
| 793 |
call back2grd (deltgather(1,L),pblindex, dtmoist(1,1,L) ,im*jm) |
| 794 |
call back2grd (delqgather(1,L),pblindex, dqmoist(1,1,L,1),im*jm) |
| 795 |
call back2grd ( cpen(1,1,L),pblindex, cpen(1,1,L) ,im*jm) |
| 796 |
call back2grd ( cldls(1,1,L),pblindex, cldls(1,1,L) ,im*jm) |
| 797 |
call back2grd (cldwater(1,1,L),pblindex,cldwater(1,1,L) ,im*jm) |
| 798 |
call back2grd ( pkzgather(1,L),pblindex, pkzgather(1,L) ,im*jm) |
| 799 |
|
| 800 |
C Diagnostics: |
| 801 |
c ------------ |
| 802 |
if(icldmas.gt.0)call back2grd(tmpgather(1,L),pblindex, |
| 803 |
. tmpgather(1,L),im*jm) |
| 804 |
if(idtrain.gt.0)call back2grd(pkegather(1,L),pblindex, |
| 805 |
. pkegather(1,L),im*jm) |
| 806 |
if(idtls.gt.0)call back2grd(deltrnev(1,L),pblindex, |
| 807 |
. deltrnev(1,L),im*jm) |
| 808 |
if(idqls.gt.0)call back2grd(delqrnev(1,L),pblindex, |
| 809 |
. delqrnev(1,L),im*jm) |
| 810 |
if(icldnp.gt.0)call back2grd(cldsr(1,1,L),pblindex, |
| 811 |
. cldsr(1,1,L),im*jm) |
| 812 |
enddo |
| 813 |
|
| 814 |
c Tracers |
| 815 |
c ------- |
| 816 |
do nt = 1,ntracer-ptracer |
| 817 |
do L = 1,lm |
| 818 |
call back2grd (delugather(1,L,nt),pblindex, |
| 819 |
. dqmoist(1,1,L,ptracer+nt),im*jm) |
| 820 |
enddo |
| 821 |
enddo |
| 822 |
|
| 823 |
|
| 824 |
C ********************************************************************** |
| 825 |
C BUMP DIAGNOSTICS |
| 826 |
C ********************************************************************** |
| 827 |
|
| 828 |
c Determine Level Indices for Low-Mid-High Cloud Regions |
| 829 |
c ------------------------------------------------------ |
| 830 |
lowlevel = lm |
| 831 |
midlevel = lm |
| 832 |
do L = lm-1,1,-1 |
| 833 |
pcheck = (1000.-ptop)*sig(l) + ptop |
| 834 |
if (pcheck.gt.700.0) lowlevel = L |
| 835 |
if (pcheck.gt.400.0) midlevel = L |
| 836 |
enddo |
| 837 |
|
| 838 |
|
| 839 |
c Clear-Sky (Above 400mb) Temperature |
| 840 |
c ----------------------------------- |
| 841 |
if( itmpuclr.ne.0 .or. isphuclr.ne.0 ) then |
| 842 |
do j = 1,jm |
| 843 |
do i = 1,im |
| 844 |
totcld(i,j) = 0.0 |
| 845 |
enddo |
| 846 |
enddo |
| 847 |
do L = 1,midlevel |
| 848 |
do j = 1,jm |
| 849 |
do i = 1,im |
| 850 |
if(cldls(i,j,L).ne.0.0.or.cpen(i,j,L).ne.0.0)totcld(i,j) = 1.0 |
| 851 |
enddo |
| 852 |
enddo |
| 853 |
enddo |
| 854 |
do L = 1,lm |
| 855 |
if( itmpuclr.ne.0 ) then |
| 856 |
do i = 1,im*jm |
| 857 |
if( totcld(i,1).eq.0.0 ) then |
| 858 |
qdiag(i,1,itmpuclr +L-1) = qdiag(i,1,itmpuclr +L-1) + |
| 859 |
. tz(i,1,L)*pkzgather(i,L) |
| 860 |
qdiag(i,1,itmpuclrc+L-1) = qdiag(i,1,itmpuclrc+L-1) + 1.0 |
| 861 |
endif |
| 862 |
enddo |
| 863 |
endif |
| 864 |
|
| 865 |
if( isphuclr.ne.0 ) then |
| 866 |
do i = 1,im*jm |
| 867 |
if( totcld(i,1).eq.0.0 ) then |
| 868 |
qdiag(i,1,isphuclr +L-1) = qdiag(i,1,isphuclr +L-1) + |
| 869 |
. qz(i,1,L,1)*1000.0 |
| 870 |
qdiag(i,1,isphuclrc+L-1) = qdiag(i,1,isphuclrc+L-1) + 1.0 |
| 871 |
endif |
| 872 |
enddo |
| 873 |
endif |
| 874 |
enddo |
| 875 |
endif |
| 876 |
|
| 877 |
c Sub-Cloud Layer |
| 878 |
c ------------------------- |
| 879 |
if( ipsubcld.ne.0 ) then |
| 880 |
do j = 1,jm |
| 881 |
do i = 1,im |
| 882 |
qdiag(i,j,ipsubcld ) = qdiag(i,j,ipsubcld ) + psubcldg (i,j) |
| 883 |
qdiag(i,j,ipsubcldc) = qdiag(i,j,ipsubcldc) + psubcldgc(i,j) |
| 884 |
enddo |
| 885 |
enddo |
| 886 |
endif |
| 887 |
|
| 888 |
c Non-Precipitating Cloud Fraction |
| 889 |
c -------------------------------- |
| 890 |
if( icldnp.ne.0 ) then |
| 891 |
do L = 1,lm |
| 892 |
do j = 1,jm |
| 893 |
do i = 1,im |
| 894 |
qdiag(i,j,icldnp+L-1) = qdiag(i,j,icldnp+L-1) + cldsr(i,j,L) |
| 895 |
enddo |
| 896 |
enddo |
| 897 |
enddo |
| 898 |
ncldnp = ncldnp + 1 |
| 899 |
endif |
| 900 |
|
| 901 |
c Moist Processes Heating Rate |
| 902 |
c ---------------------------- |
| 903 |
if(imoistt.gt.0) then |
| 904 |
do L = 1,lm |
| 905 |
do i = 1,im*jm |
| 906 |
qdiag(i,1,imoistt+L-1) = qdiag(i,1,imoistt+L-1) + |
| 907 |
. (dtmoist(i,1,L)*sday*pkzgather(i,L)/pz(i,1)) |
| 908 |
enddo |
| 909 |
enddo |
| 910 |
endif |
| 911 |
|
| 912 |
c Moist Processes Moistening Rate |
| 913 |
c ------------------------------- |
| 914 |
if(imoistq.gt.0) then |
| 915 |
do L = 1,lm |
| 916 |
do j = 1,jm |
| 917 |
do i = 1,im |
| 918 |
qdiag(i,j,imoistq+L-1) = qdiag(i,j,imoistq+L-1) + |
| 919 |
. (dqmoist(i,j,L,1)*sday*1000.0/pz(i,j)) |
| 920 |
enddo |
| 921 |
enddo |
| 922 |
enddo |
| 923 |
endif |
| 924 |
|
| 925 |
c Cloud Mass Flux |
| 926 |
c --------------- |
| 927 |
if(icldmas.gt.0) then |
| 928 |
do L = 1,lm |
| 929 |
do i = 1,im*jm |
| 930 |
qdiag(i,1,icldmas+L-1) = qdiag(i,1,icldmas+L-1) + tmpgather(i,L) |
| 931 |
enddo |
| 932 |
enddo |
| 933 |
endif |
| 934 |
|
| 935 |
c Detrained Cloud Mass Flux |
| 936 |
c ------------------------- |
| 937 |
if(idtrain.gt.0) then |
| 938 |
do L = 1,lm |
| 939 |
do i = 1,im*jm |
| 940 |
qdiag(i,1,idtrain+L-1) = qdiag(i,1,idtrain+L-1) + pkegather(i,L) |
| 941 |
enddo |
| 942 |
enddo |
| 943 |
endif |
| 944 |
|
| 945 |
c Grid-Scale Condensational Heating Rate |
| 946 |
c -------------------------------------- |
| 947 |
if(idtls.gt.0) then |
| 948 |
do L = 1,lm |
| 949 |
do i = 1,im*jm |
| 950 |
qdiag(i,1,idtls+L-1) = qdiag(i,1,idtls+L-1) + deltrnev(i,L) |
| 951 |
enddo |
| 952 |
enddo |
| 953 |
endif |
| 954 |
|
| 955 |
c Grid-Scale Condensational Moistening Rate |
| 956 |
c ----------------------------------------- |
| 957 |
if(idqls.gt.0) then |
| 958 |
do L = 1,lm |
| 959 |
do i = 1,im*jm |
| 960 |
qdiag(i,1,idqls+L-1) = qdiag(i,1,idqls+L-1) + delqrnev(i,L) |
| 961 |
enddo |
| 962 |
enddo |
| 963 |
endif |
| 964 |
|
| 965 |
c Total Precipitation |
| 966 |
c ------------------- |
| 967 |
if(ipreacc.gt.0) then |
| 968 |
do j = 1,jm |
| 969 |
do i = 1,im |
| 970 |
qdiag(i,j,ipreacc) = qdiag(i,j,ipreacc) |
| 971 |
. + ( lsp_new(I,j) |
| 972 |
. + snow_new(I,j) |
| 973 |
. + conv_new(i,j) ) *sday*tminv |
| 974 |
enddo |
| 975 |
enddo |
| 976 |
endif |
| 977 |
|
| 978 |
c Convective Precipitation |
| 979 |
c ------------------------ |
| 980 |
if(iprecon.gt.0) then |
| 981 |
do i = 1,im*jm |
| 982 |
qdiag(i,1,iprecon) = qdiag(i,1,iprecon) + raincgath(i)*sday*tminv |
| 983 |
enddo |
| 984 |
endif |
| 985 |
|
| 986 |
C ********************************************************************** |
| 987 |
C **** Fill Rainfall and Snowfall Arrays for Land Surface Model **** |
| 988 |
C **** Note: Precip Rates work when DT(turb)<DT(moist) **** |
| 989 |
C ********************************************************************** |
| 990 |
|
| 991 |
do j = 1,jm |
| 992 |
do i = 1,im |
| 993 |
rainlsp (i,j) = rainlsp (i,j) + lsp_new(i,j)*tminv |
| 994 |
rainconv(i,j) = rainconv(i,j) + conv_new(i,j)*tminv |
| 995 |
snowfall(i,j) = snowfall(i,j) + snow_new(i,j)*tminv |
| 996 |
enddo |
| 997 |
enddo |
| 998 |
|
| 999 |
C ********************************************************************** |
| 1000 |
C *** Compute Time-averaged Quantities for Radiation *** |
| 1001 |
C *** CPEN => Cloud Fraction from RAS *** |
| 1002 |
C *** CLDLS => Cloud Fraction from RNEVP *** |
| 1003 |
C ********************************************************************** |
| 1004 |
|
| 1005 |
do j = 1,jm |
| 1006 |
do i = 1,im |
| 1007 |
cldhi (i,j) = 0. |
| 1008 |
cldmid(i,j) = 0. |
| 1009 |
cldlow(i,j) = 0. |
| 1010 |
cldtmp(i,j) = 0. |
| 1011 |
cldprs(i,j) = 0. |
| 1012 |
tmpimjm(i,j) = 0. |
| 1013 |
enddo |
| 1014 |
enddo |
| 1015 |
|
| 1016 |
c Set Moist-Process Memory Coefficient |
| 1017 |
c ------------------------------------ |
| 1018 |
cldras_mem = 1.0-tmstp/ 3600.0 |
| 1019 |
cldlsp_mem = 1.0-tmstp/(3600.0*3) |
| 1020 |
|
| 1021 |
do L = 1,lm |
| 1022 |
do i = 1,im*jm |
| 1023 |
plev = sig(L)*pz(i,1)+ptop |
| 1024 |
|
| 1025 |
c Compute Time-averaged Cloud and Water Amounts for Longwave Radiation |
| 1026 |
c -------------------------------------------------------------------- |
| 1027 |
watnow = cldwater(i,1,L) |
| 1028 |
if( plev.le.500.0 ) then |
| 1029 |
cldras = min( max( cldras_lw(i,1,L)*cldras_mem,cpen(i,1,L)),1.0) |
| 1030 |
else |
| 1031 |
cldras = 0.0 |
| 1032 |
endif |
| 1033 |
cldlsp = min( max( cldlsp_lw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0) |
| 1034 |
|
| 1035 |
if( cldras.lt.cldmin ) cldras = 0.0 |
| 1036 |
if( cldlsp.lt.cldmin ) cldlsp = 0.0 |
| 1037 |
|
| 1038 |
cldnow = max( cldlsp,cldras ) |
| 1039 |
|
| 1040 |
lwlz(i,1,L) = ( nlwlz*lwlz(i,1,L) + watnow)/(nlwlz +1) |
| 1041 |
cldtot_lw(i,1,L) = (nlwcld*cldtot_lw(i,1,L) + cldnow)/(nlwcld+1) |
| 1042 |
cldlsp_lw(i,1,L) = (nlwcld*cldlsp_lw(i,1,L) + cldlsp)/(nlwcld+1) |
| 1043 |
cldras_lw(i,1,L) = (nlwcld*cldras_lw(i,1,L) + cldras)/(nlwcld+1) |
| 1044 |
|
| 1045 |
|
| 1046 |
c Compute Time-averaged Cloud and Water Amounts for Shortwave Radiation |
| 1047 |
c --------------------------------------------------------------------- |
| 1048 |
watnow = cldwater(i,1,L) |
| 1049 |
if( plev.le.500.0 ) then |
| 1050 |
cldras = min( max(cldras_sw(i,1,L)*cldras_mem, cpen(i,1,L)),1.0) |
| 1051 |
else |
| 1052 |
cldras = 0.0 |
| 1053 |
endif |
| 1054 |
cldlsp = min( max(cldlsp_sw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0) |
| 1055 |
|
| 1056 |
if( cldras.lt.cldmin ) cldras = 0.0 |
| 1057 |
if( cldlsp.lt.cldmin ) cldlsp = 0.0 |
| 1058 |
|
| 1059 |
cldnow = max( cldlsp,cldras ) |
| 1060 |
|
| 1061 |
swlz(i,1,L) = ( nswlz*swlz(i,1,L) + watnow)/(nswlz +1) |
| 1062 |
cldtot_sw(i,1,L) = (nswcld*cldtot_sw(i,1,L) + cldnow)/(nswcld+1) |
| 1063 |
cldlsp_sw(i,1,L) = (nswcld*cldlsp_sw(i,1,L) + cldlsp)/(nswcld+1) |
| 1064 |
cldras_sw(i,1,L) = (nswcld*cldras_sw(i,1,L) + cldras)/(nswcld+1) |
| 1065 |
|
| 1066 |
|
| 1067 |
c Compute Instantaneous Low-Mid-High Maximum Overlap Cloud Fractions |
| 1068 |
c ---------------------------------------------------------------------- |
| 1069 |
|
| 1070 |
if( L.lt.midlevel ) cldhi (i,1) = max( cldnow,cldhi (i,1) ) |
| 1071 |
if( L.ge.midlevel .and. |
| 1072 |
. L.lt.lowlevel ) cldmid(i,1) = max( cldnow,cldmid(i,1) ) |
| 1073 |
if( L.ge.lowlevel ) cldlow(i,1) = max( cldnow,cldlow(i,1) ) |
| 1074 |
|
| 1075 |
c Compute Cloud-Top Temperature and Pressure |
| 1076 |
c ------------------------------------------ |
| 1077 |
cldtmp(i,1) = cldtmp(i,1) + cldnow*pkzgather(i,L) |
| 1078 |
. * ( tz(i,1,L) + dtmoist(i,1,L)*tmstp/pz(i,1) ) |
| 1079 |
cldprs(i,1) = cldprs(i,1) + cldnow*plev |
| 1080 |
tmpimjm(i,1) = tmpimjm(i,1) + cldnow |
| 1081 |
|
| 1082 |
enddo |
| 1083 |
enddo |
| 1084 |
|
| 1085 |
c Compute Instantanious Total 2-D Cloud Fraction |
| 1086 |
c ---------------------------------------------- |
| 1087 |
do j = 1,jm |
| 1088 |
do i = 1,im |
| 1089 |
totcld(i,j) = 1.0 - (1.-cldhi (i,j)) |
| 1090 |
. * (1.-cldmid(i,j)) |
| 1091 |
. * (1.-cldlow(i,j)) |
| 1092 |
enddo |
| 1093 |
enddo |
| 1094 |
|
| 1095 |
|
| 1096 |
C ********************************************************************** |
| 1097 |
C *** Fill Cloud Top Pressure and Temperature Diagnostic *** |
| 1098 |
C ********************************************************************** |
| 1099 |
|
| 1100 |
if(icldtmp.gt.0) then |
| 1101 |
do j = 1,jm |
| 1102 |
do i = 1,im |
| 1103 |
if( cldtmp(i,j).gt.0.0 ) then |
| 1104 |
qdiag(i,j,icldtmp) = qdiag(i,j,icldtmp) + |
| 1105 |
. cldtmp(i,j)*totcld(i,j)/tmpimjm(i,j) |
| 1106 |
qdiag(i,j,icttcnt) = qdiag(i,j,icttcnt) + totcld(i,j) |
| 1107 |
endif |
| 1108 |
enddo |
| 1109 |
enddo |
| 1110 |
endif |
| 1111 |
|
| 1112 |
if(icldprs.gt.0) then |
| 1113 |
do j = 1,jm |
| 1114 |
do i = 1,im |
| 1115 |
if( cldprs(i,j).gt.0.0 ) then |
| 1116 |
qdiag(i,j,icldprs) = qdiag(i,j,icldprs) + |
| 1117 |
. cldprs(i,j)*totcld(i,j)/tmpimjm(i,j) |
| 1118 |
qdiag(i,j,ictpcnt) = qdiag(i,j,ictpcnt) + totcld(i,j) |
| 1119 |
endif |
| 1120 |
enddo |
| 1121 |
enddo |
| 1122 |
endif |
| 1123 |
|
| 1124 |
C ********************************************************************** |
| 1125 |
C **** INCREMENT COUNTERS **** |
| 1126 |
C ********************************************************************** |
| 1127 |
|
| 1128 |
nlwlz = nlwlz + 1 |
| 1129 |
nswlz = nswlz + 1 |
| 1130 |
|
| 1131 |
nlwcld = nlwcld + 1 |
| 1132 |
nswcld = nswcld + 1 |
| 1133 |
|
| 1134 |
nmoistt = nmoistt + 1 |
| 1135 |
nmoistq = nmoistq + 1 |
| 1136 |
npreacc = npreacc + 1 |
| 1137 |
nprecon = nprecon + 1 |
| 1138 |
|
| 1139 |
ncldmas = ncldmas + 1 |
| 1140 |
ndtrain = ndtrain + 1 |
| 1141 |
|
| 1142 |
ndtls = ndtls + 1 |
| 1143 |
ndqls = ndqls + 1 |
| 1144 |
|
| 1145 |
RETURN |
| 1146 |
END |
| 1147 |
SUBROUTINE RAS( NN, LEN, LENC, K, NLTOP, nlayr, DT |
| 1148 |
*, UOI, ntracer, POI, QOI, PRS, PRJ, rnd, ncrnd |
| 1149 |
*, RAINS, CLN, CLF, cldmas, detrain |
| 1150 |
*, cp,grav,rkappa,alhl,rhfrac,rasmax ) |
| 1151 |
C |
| 1152 |
C********************************************************************* |
| 1153 |
C*********************** ARIES MODEL ******************************* |
| 1154 |
C********************* SUBROUTINE RAS ***************************** |
| 1155 |
C********************** 16 MARCH 1988 ****************************** |
| 1156 |
C********************************************************************* |
| 1157 |
C |
| 1158 |
PARAMETER (KRMIN=01) |
| 1159 |
PARAMETER (ICM=1000) |
| 1160 |
PARAMETER (CMB2PA=100.0) |
| 1161 |
PARAMETER (rknob = 10.) |
| 1162 |
C |
| 1163 |
integer ntracer |
| 1164 |
integer nltop,nlayr |
| 1165 |
DIMENSION UOI(len,nlayr,ntracer), POI(len,K) |
| 1166 |
DIMENSION QOI(len,K), PRS(len,K+1), PRJ(len,K+1) |
| 1167 |
dimension rnd(ncrnd) |
| 1168 |
C |
| 1169 |
DIMENSION RAINS(len,K), CLN(len,K), CLF(len,K) |
| 1170 |
DIMENSION cldmas(len,K), detrain(len,K) |
| 1171 |
DIMENSION TCU(len,K), QCU(len,K) |
| 1172 |
real ucu(len,K,ntracer) |
| 1173 |
DIMENSION ALF(len,K), BET(len,K), GAM(len,K) |
| 1174 |
*, ETA(len,K), HOI(len,K) |
| 1175 |
*, PRH(len,K), PRI(len,K) |
| 1176 |
DIMENSION HST(len,K), QOL(len,K), GMH(len,K) |
| 1177 |
|
| 1178 |
DIMENSION TX1(len), TX2(len), TX3(len), TX4(len), TX5(len) |
| 1179 |
*, TX6(len), TX7(len), TX8(len), TX9(len) |
| 1180 |
*, TX11(len), TX12(len), TX13(len), TX14(len,ntracer) |
| 1181 |
*, TX15(len), TX16(len) |
| 1182 |
*, WFN(len), IA1(len), IA2(len), IA3(len) |
| 1183 |
DIMENSION cloudn(len), pcu(len) |
| 1184 |
|
| 1185 |
real rhfrac(len),rasmax |
| 1186 |
|
| 1187 |
DIMENSION IC(ICM), IRND(icm) |
| 1188 |
dimension cmass(len,K) |
| 1189 |
LOGICAL SETRAS |
| 1190 |
|
| 1191 |
do L = 1,k |
| 1192 |
do I = 1,LENC |
| 1193 |
rains(i,l) = 0. |
| 1194 |
enddo |
| 1195 |
enddo |
| 1196 |
|
| 1197 |
p00 = 1000. |
| 1198 |
crtmsf = 0. |
| 1199 |
|
| 1200 |
C The numerator here is the fraction of the subcloud layer mass flux |
| 1201 |
C allowed to entrain into the cloud |
| 1202 |
|
| 1203 |
CCC FRAC = 1./dt |
| 1204 |
FRAC = 0.5/dt |
| 1205 |
|
| 1206 |
KM1 = K - 1 |
| 1207 |
KP1 = K + 1 |
| 1208 |
C we want the ras adjustment time scale to be one hour (indep of dt) |
| 1209 |
RASBLF = 1./3600. |
| 1210 |
C |
| 1211 |
KPRV = KM1 |
| 1212 |
C Removed KRMAX parameter |
| 1213 |
KCR = MIN(KM1,nlayr-2) |
| 1214 |
KFX = KM1 - KCR |
| 1215 |
NCMX = KFX + NCRND |
| 1216 |
C |
| 1217 |
IF (KFX .GT. 0) THEN |
| 1218 |
DO NC=1,KFX |
| 1219 |
IC(NC) = K - NC |
| 1220 |
ENDDO |
| 1221 |
ENDIF |
| 1222 |
C |
| 1223 |
IF (NCRND .GT. 0) THEN |
| 1224 |
DO I=1,ncrnd |
| 1225 |
IRND(I) = (RND(I)-0.0005)*(KCR-KRMIN+1) |
| 1226 |
IRND(I) = IRND(I) + KRMIN |
| 1227 |
ENDDO |
| 1228 |
C |
| 1229 |
DO NC=1,NCRND |
| 1230 |
IC(KFX+NC) = IRND(NC) |
| 1231 |
ENDDO |
| 1232 |
ENDIF |
| 1233 |
C |
| 1234 |
DO 100 NC=1,NCMX |
| 1235 |
C |
| 1236 |
IF (NC .EQ. 1 ) THEN |
| 1237 |
SETRAS = .TRUE. |
| 1238 |
ELSE |
| 1239 |
SETRAS = .FALSE. |
| 1240 |
ENDIF |
| 1241 |
IB = IC(NC) |
| 1242 |
|
| 1243 |
c Initialize Cloud Fraction Array |
| 1244 |
c ------------------------------- |
| 1245 |
do i = 1,lenc |
| 1246 |
cloudn(i) = 0.0 |
| 1247 |
enddo |
| 1248 |
|
| 1249 |
CALL CLOUD(nn,LEN, LENC, K, NLTOP, nlayr, IB, RASBLF,SETRAS,FRAC |
| 1250 |
*, CP, ALHL, RKAPPA, GRAV, P00, CRTMSF |
| 1251 |
*, POI, QOI, UOI, Ntracer, PRS, PRJ |
| 1252 |
*, PCU, CLOUDN, TCU, QCU, UCU, CMASS |
| 1253 |
*, ALF, BET, GAM, PRH, PRI, HOI, ETA |
| 1254 |
*, HST, QOL, GMH |
| 1255 |
*, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9 |
| 1256 |
*, WFN, TX11, TX12, TX13, TX14, TX15 |
| 1257 |
*, IA1,IA2,IA3,rhfrac) |
| 1258 |
|
| 1259 |
C Compute fraction of grid box into which rain re-evap occurs (clf) |
| 1260 |
c ----------------------------------------------------------------- |
| 1261 |
do i = 1,lenc |
| 1262 |
|
| 1263 |
c mass in detrainment layer |
| 1264 |
c ------------------------- |
| 1265 |
tx1(i) = cmb2pa * (prs(i,ib+1) - prs(i,ib))/(grav*dt) |
| 1266 |
|
| 1267 |
c ratio of detraining cloud mass to mass in detrainment layer |
| 1268 |
c ----------------------------------------------------------- |
| 1269 |
tx1(i) = rhfrac(i)*rknob * cmass(i,ib) / tx1(i) |
| 1270 |
if(cmass(i,K).gt.0.) clf(i,ib) = clf(i,ib) + tx1(i) |
| 1271 |
if( clf(i,ib).gt.1.) clf(i,ib) = 1. |
| 1272 |
enddo |
| 1273 |
|
| 1274 |
c Compute Total Cloud Mass Flux |
| 1275 |
c ***************************** |
| 1276 |
do L=ib,k |
| 1277 |
do i=1,lenc |
| 1278 |
cmass(i,L) = rhfrac(i)*cmass(i,L) * dt |
| 1279 |
enddo |
| 1280 |
enddo |
| 1281 |
|
| 1282 |
do L=ib,k |
| 1283 |
do i=1,lenc |
| 1284 |
cldmas(i,L) = cldmas(i,L) + cmass(i,L) |
| 1285 |
enddo |
| 1286 |
enddo |
| 1287 |
|
| 1288 |
do i=1,lenc |
| 1289 |
detrain(i,ib) = detrain(i,ib) + cmass(i,ib) |
| 1290 |
enddo |
| 1291 |
|
| 1292 |
DO L=IB,K |
| 1293 |
DO I=1,LENC |
| 1294 |
POI(I,L) = POI(I,L) + TCU(I,L) * DT * rhfrac(i) |
| 1295 |
QOI(I,L) = QOI(I,L) + QCU(I,L) * DT * rhfrac(i) |
| 1296 |
ENDDO |
| 1297 |
ENDDO |
| 1298 |
DO NT=1,Ntracer |
| 1299 |
DO L=IB,K |
| 1300 |
DO I=1,LENC |
| 1301 |
UOI(I,L+nltop-1,NT)=UOI(I,L+nltop-1,NT)+UCU(I,L,NT)*DT*rhfrac(i) |
| 1302 |
ENDDO |
| 1303 |
ENDDO |
| 1304 |
ENDDO |
| 1305 |
DO I=1,LENC |
| 1306 |
rains(I,ib) = rains(I,ib) + PCU(I)*dt * rhfrac(i) |
| 1307 |
ENDDO |
| 1308 |
|
| 1309 |
100 CONTINUE |
| 1310 |
|
| 1311 |
c Fill Convective Cloud Fractions based on 3-D Rain Amounts |
| 1312 |
c --------------------------------------------------------- |
| 1313 |
do L=k-1,1,-1 |
| 1314 |
do i=1,lenc |
| 1315 |
tx1(i) = 100*(prs(i,L+1)-prs(i,L))/grav |
| 1316 |
cln(i,L) = min(1600*rains(i,L)/tx1(i),rasmax ) |
| 1317 |
enddo |
| 1318 |
enddo |
| 1319 |
|
| 1320 |
RETURN |
| 1321 |
END |
| 1322 |
|
| 1323 |
subroutine rndcloud (iras,nrnd,rnd,myid) |
| 1324 |
implicit none |
| 1325 |
integer n,iras,nrnd,myid |
| 1326 |
real random_numbx |
| 1327 |
real rnd(nrnd) |
| 1328 |
integer irm |
| 1329 |
parameter (irm = 1000) |
| 1330 |
real random(irm) |
| 1331 |
integer i,mcheck,numrand,iseed,index |
| 1332 |
logical first |
| 1333 |
data first /.true./ |
| 1334 |
integer iras0 |
| 1335 |
data iras0 /0/ |
| 1336 |
save random, iras0 |
| 1337 |
|
| 1338 |
if(nrnd.eq.0.)then |
| 1339 |
do i = 1,nrnd |
| 1340 |
rnd(i) = 0 |
| 1341 |
enddo |
| 1342 |
if(first .and. myid.eq.0) print *,' NO RANDOM CLOUDS IN RAS ' |
| 1343 |
go to 100 |
| 1344 |
endif |
| 1345 |
|
| 1346 |
mcheck = mod(iras-1,irm/nrnd) |
| 1347 |
|
| 1348 |
c First Time In From a Continuing RESTART (IRAS.GT.1) or Reading a New RESTART |
| 1349 |
c ---------------------------------------------------------------------------- |
| 1350 |
if( first.and.(iras.gt.1) .or. iras.ne.iras0+1 )then |
| 1351 |
if( myid.eq.0 ) print *, 'Recreating Rand Numb Array in RNDCLOUD' |
| 1352 |
if( myid.eq.0 ) print *, 'IRAS: ',iras,' IRAS0: ',iras0 |
| 1353 |
numrand = mod(iras,irm/nrnd) * nrnd |
| 1354 |
iseed = iras * nrnd - numrand |
| 1355 |
call random_seedx(iseed) |
| 1356 |
do i = 1,irm |
| 1357 |
random(i) = random_numbx() |
| 1358 |
enddo |
| 1359 |
index = (iras-1)*nrnd |
| 1360 |
|
| 1361 |
c Multiple Time In But have Used Up all 1000 numbers (MCHECK.EQ.0) |
| 1362 |
c ---------------------------------------------------------------- |
| 1363 |
else if (mcheck.eq.0) then |
| 1364 |
iseed = (iras-1)*nrnd |
| 1365 |
call random_seedx(iseed) |
| 1366 |
do i = 1,irm |
| 1367 |
random(i) = random_numbx() |
| 1368 |
enddo |
| 1369 |
index = iseed |
| 1370 |
|
| 1371 |
c Multiple Time In But have NOT Used Up all 1000 numbers (MCHECK.NE.0) |
| 1372 |
c -------------------------------------------------------------------- |
| 1373 |
else |
| 1374 |
index = (iras-1)*nrnd |
| 1375 |
endif |
| 1376 |
|
| 1377 |
index = mod(index,irm) |
| 1378 |
if( index+nrnd.gt.1000 ) index=1000-nrnd |
| 1379 |
|
| 1380 |
do n = 1,nrnd |
| 1381 |
rnd(n) = random(index+n) |
| 1382 |
enddo |
| 1383 |
|
| 1384 |
100 continue |
| 1385 |
first = .false. |
| 1386 |
iras0 = iras |
| 1387 |
return |
| 1388 |
end |
| 1389 |
|
| 1390 |
real function random_numbx() |
| 1391 |
implicit none |
| 1392 |
#if CRAY |
| 1393 |
real ranf |
| 1394 |
random_numbx = ranf() |
| 1395 |
#endif |
| 1396 |
#if SGI |
| 1397 |
real rand |
| 1398 |
random_numbx = rand() |
| 1399 |
#endif |
| 1400 |
return |
| 1401 |
end |
| 1402 |
subroutine random_seedx (iseed) |
| 1403 |
implicit none |
| 1404 |
integer iseed |
| 1405 |
#if CRAY |
| 1406 |
call ranset (iseed) |
| 1407 |
#endif |
| 1408 |
#if SGI |
| 1409 |
integer*4 seed |
| 1410 |
seed = iseed |
| 1411 |
call srand (seed) |
| 1412 |
#endif |
| 1413 |
return |
| 1414 |
end |
| 1415 |
|
| 1416 |
SUBROUTINE CLOUD(nn,LEN, LENC, K, NLTOP, nlayr, IC, RASALF, |
| 1417 |
*, SETRAS, FRAC |
| 1418 |
*, CP, ALHL, RKAP, GRAV, P00, CRTMSF |
| 1419 |
*, POI, QOI, UOI, Ntracer, PRS, PRJ |
| 1420 |
*, PCU, CLN, TCU, QCU, UCU, CMASS |
| 1421 |
*, ALF, BET, GAM, PRH, PRI, HOL, ETA |
| 1422 |
*, HST, QOL, GMH |
| 1423 |
*, TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, ALM |
| 1424 |
*, WFN, AKM, QS1, CLF, UHT, WLQ |
| 1425 |
*, IA, I1, I2,rhfrac) |
| 1426 |
C |
| 1427 |
C********************************************************************* |
| 1428 |
C******************** Relaxed Arakawa-Schubert *********************** |
| 1429 |
C********************* Plug Compatible Version ********************** |
| 1430 |
C************************ SUBROUTINE CLOUD *************************** |
| 1431 |
C************************* 23 JULY 1992 *************************** |
| 1432 |
C********************************************************************* |
| 1433 |
C********************************************************************* |
| 1434 |
C********************************************************************* |
| 1435 |
C************************** Developed By ***************************** |
| 1436 |
C************************** ***************************** |
| 1437 |
C************************ Shrinivas Moorthi ************************** |
| 1438 |
C************************ and ************************** |
| 1439 |
C************************ Max J. Suarez ***************************** |
| 1440 |
C************************ ***************************** |
| 1441 |
C******************** Laboratory for Atmospheres ********************* |
| 1442 |
C****************** NASA/GSFC, Greenbelt, MD 20771 ******************* |
| 1443 |
C********************************************************************* |
| 1444 |
C********************************************************************* |
| 1445 |
C |
| 1446 |
C The calculations of Moorthi and Suarez (1992, MWR) are |
| 1447 |
C contained in the CLOUD routine. |
| 1448 |
C It is probably advisable, at least initially, to treat CLOUD |
| 1449 |
C as a black box that computes the single cloud adjustments. RAS, |
| 1450 |
C on the other hand, can be tailored to each GCMs configuration |
| 1451 |
C (ie, number and placement of levels, nature of boundary layer, |
| 1452 |
C time step and frequency with which RAS is called). |
| 1453 |
C |
| 1454 |
C |
| 1455 |
C Input: |
| 1456 |
C ------ |
| 1457 |
C |
| 1458 |
C LEN : The inner dimension of update and input arrays. |
| 1459 |
C |
| 1460 |
C LENC : The run: the number of soundings processes in a single call. |
| 1461 |
C RAS works on the first LENC of the LEN soundings |
| 1462 |
C passed. This allows working on pieces of the world |
| 1463 |
C say for multitasking, without declaring temporary arrays |
| 1464 |
C and copying the data to and from them. This is an f77 |
| 1465 |
C version. An F90 version would have to allow more |
| 1466 |
C flexibility in the argument declarations. Obviously |
| 1467 |
C (LENC<=LEN). |
| 1468 |
C |
| 1469 |
C K : Number of vertical layers (increasing downwards). |
| 1470 |
C Need not be the same as the number of layers in the |
| 1471 |
C GCM, since it is the outer dimension. The bottom layer |
| 1472 |
C (K) is the subcloud layer. |
| 1473 |
C |
| 1474 |
C IC : Detrainment level to check for presence of convection |
| 1475 |
C |
| 1476 |
C RASALF : Relaxation parameter (< 1.) for present cloud-type |
| 1477 |
C |
| 1478 |
C SETRAS : Logical parameter to control re-calculation of |
| 1479 |
C saturation specific humidity and mid level P**kappa |
| 1480 |
C |
| 1481 |
C FRAC : Fraction of the PBL (layer K) mass allowed to be used |
| 1482 |
C by a cloud-type in time DT |
| 1483 |
C |
| 1484 |
C CP : Specific heat at constant pressure |
| 1485 |
C |
| 1486 |
C ALHL : Latent Heat of condensation |
| 1487 |
C |
| 1488 |
C RKAP : R/Cp, where R is the gas constant |
| 1489 |
C |
| 1490 |
C GRAV : Acceleration due to gravity |
| 1491 |
C |
| 1492 |
C P00 : A reference pressure in hPa, useually 1000 hPa |
| 1493 |
C |
| 1494 |
C CRTMSF : Critical value of mass flux above which cloudiness at |
| 1495 |
C the detrainment layer of that cloud-type is assumed. |
| 1496 |
C Affects only cloudiness calculation. |
| 1497 |
C |
| 1498 |
C POI : 2D array of dimension (LEN,K) containing potential |
| 1499 |
C temperature. Updated but not initialized by RAS. |
| 1500 |
C |
| 1501 |
C QOI : 2D array of dimension (LEN,K) containing specific |
| 1502 |
C humidity. Updated but not initialized by RAS. |
| 1503 |
C |
| 1504 |
C UOI : 3D array of dimension (LEN,K,NTRACER) containing tracers |
| 1505 |
C Updated but not initialized by RAS. |
| 1506 |
C |
| 1507 |
C PRS : 2D array of dimension (LEN,K+1) containing pressure |
| 1508 |
C in hPa at the interfaces of K-layers from top of the |
| 1509 |
C atmosphere to the bottom. Not modified. |
| 1510 |
C |
| 1511 |
C PRJ : 2D array of dimension (LEN,K+1) containing (PRS/P00) ** |
| 1512 |
C RKAP. i.e. Exner function at layer edges. Not modified. |
| 1513 |
C |
| 1514 |
C rhfrac : 1D array of dimension (LEN) containing a rel.hum. scaling |
| 1515 |
C fraction. Not modified. |
| 1516 |
C |
| 1517 |
C Output: |
| 1518 |
C ------- |
| 1519 |
C |
| 1520 |
C PCU : 1D array of length LEN containing accumulated |
| 1521 |
C precipitation in mm/sec. |
| 1522 |
C |
| 1523 |
C CLN : 2D array of dimension (LEN,K) containing cloudiness |
| 1524 |
C Note: CLN is bumped but NOT initialized |
| 1525 |
C |
| 1526 |
C TCU : 2D array of dimension (LEN,K) containing accumulated |
| 1527 |
C convective heating (K/sec). |
| 1528 |
C |
| 1529 |
C QCU : 2D array of dimension (LEN,K) containing accumulated |
| 1530 |
C convective drying (kg/kg/sec). |
| 1531 |
C |
| 1532 |
C CMASS : 2D array of dimension (LEN,K) containing the |
| 1533 |
C cloud mass flux (kg/sec). Filled from cloud top |
| 1534 |
C to base. |
| 1535 |
C |
| 1536 |
C Temporaries: |
| 1537 |
C |
| 1538 |
C ALF, BET, GAM, ETA, PRH, PRI, HOI, HST, QOL, GMH are temporary |
| 1539 |
C 2D real arrays of dimension of at least (LENC,K) where LENC is |
| 1540 |
C the horizontal dimension over which convection is invoked. |
| 1541 |
C |
| 1542 |
C |
| 1543 |
C TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9, AKM, QS1, CLF, UHT |
| 1544 |
C VHT, WLQ WFN are temporary real arrays of length at least LENC |
| 1545 |
C |
| 1546 |
C IA, I1, and I2 are temporary integer arrays of length LENC |
| 1547 |
C |
| 1548 |
C |
| 1549 |
C************************************************************************ |
| 1550 |
C |
| 1551 |
C |
| 1552 |
|
| 1553 |
PARAMETER (DAYLEN=86400.0, HALF=0.5, ONE=1.0, ZERO=0.0) |
| 1554 |
PARAMETER (CMB2PA=100.0) |
| 1555 |
PARAMETER (RHMAX=0.9999) |
| 1556 |
C |
| 1557 |
integer nltop,ntracer,nlayr |
| 1558 |
DIMENSION POI(LEN,K), QOI(LEN,K), PRS(LEN,K+1) |
| 1559 |
*, PRJ(LEN,K+1) |
| 1560 |
*, TCU(LEN,K), QCU(LEN,K), CMASS(LEN,K), CLN(LEN) |
| 1561 |
real uoi(len,nlayr,ntracer) |
| 1562 |
DIMENSION ALF(LEN,K), BET(LEN,K), GAM(LEN,K) |
| 1563 |
*, PRH(LEN,K), PRI(LEN,K) |
| 1564 |
DIMENSION AKM(LENC), WFN(LENC) |
| 1565 |
DIMENSION HOL(LENC,K), QOL(LENC,K), ETA(LENC,K), HST(LENC,K) |
| 1566 |
*, GMH(LENC,K), ALM(LENC), WLQ(LENC), QS1(LENC) |
| 1567 |
*, TX1(LENC), TX2(LENC), TX3(LENC), TX4(LENC) |
| 1568 |
*, TX5(LENC), TX6(LENC), TX7(LENC), TX8(LENC) |
| 1569 |
*, CLF(LENC), PCU(LENC) |
| 1570 |
DIMENSION IA(LENC), I1(LENC), I2(LENC) |
| 1571 |
real rhfrac(len) |
| 1572 |
real ucu(len,k,ntracer),uht(len,ntracer) |
| 1573 |
LOGICAL SETRAS |
| 1574 |
|
| 1575 |
integer nt |
| 1576 |
|
| 1577 |
c Explicit Inline Directives |
| 1578 |
c -------------------------- |
| 1579 |
#if CRAY |
| 1580 |
#if f77 |
| 1581 |
cfpp$ expand (qsat) |
| 1582 |
#endif |
| 1583 |
#endif |
| 1584 |
|
| 1585 |
RKAPP1 = 1.0 + RKAP |
| 1586 |
ONEBCP = 1.0 / CP |
| 1587 |
ALBCP = ALHL * ONEBCP |
| 1588 |
ONEBG = 1.0 / GRAV |
| 1589 |
CPBG = CP * ONEBG |
| 1590 |
TWOBAL = 2.0 / ALHL |
| 1591 |
C |
| 1592 |
KM1 = K - 1 |
| 1593 |
IC1 = IC + 1 |
| 1594 |
C |
| 1595 |
C SETTIING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T. |
| 1596 |
C |
| 1597 |
|
| 1598 |
IF (SETRAS) THEN |
| 1599 |
|
| 1600 |
DO 2050 L=1,K |
| 1601 |
DO 2030 I=1,LENC |
| 1602 |
PRH(I,L) = (PRJ(I,L+1)*PRS(I,L+1) - PRJ(I,L)*PRS(I,L)) |
| 1603 |
* / ((PRS(I,L+1)-PRS(I,L)) * RKAPP1) |
| 1604 |
2030 CONTINUE |
| 1605 |
2050 CONTINUE |
| 1606 |
|
| 1607 |
DO 2070 L=1,K |
| 1608 |
DO 2060 I=1,LENC |
| 1609 |
TX5(I) = POI(I,L) * PRH(I,L) |
| 1610 |
TX1(I) = (PRS(I,L) + PRS(I,L+1)) * 0.5 |
| 1611 |
TX3(I) = TX5(I) |
| 1612 |
CALL QSAT(TX3(I), TX1(I), TX2(I), TX4(I), .TRUE.) |
| 1613 |
ALF(I,L) = TX2(I) - TX4(I) * TX5(I) |
| 1614 |
BET(I,L) = TX4(I) * PRH(I,L) |
| 1615 |
GAM(I,L) = 1.0 / ((1.0 + TX4(I)*ALBCP) * PRH(I,L)) |
| 1616 |
PRI(I,L) = (CP/CMB2PA) / (PRS(I,L+1) - PRS(I,L)) |
| 1617 |
2060 CONTINUE |
| 1618 |
2070 CONTINUE |
| 1619 |
|
| 1620 |
ENDIF |
| 1621 |
C |
| 1622 |
C |
| 1623 |
DO 10 L=1,K |
| 1624 |
DO 10 I=1,LEN |
| 1625 |
TCU(I,L) = 0.0 |
| 1626 |
QCU(I,L) = 0.0 |
| 1627 |
CMASS(I,L) = 0.0 |
| 1628 |
10 CONTINUE |
| 1629 |
|
| 1630 |
do nt = 1,ntracer |
| 1631 |
do L=1,K |
| 1632 |
do I=1,LENC |
| 1633 |
ucu(I,L,nt) = 0.0 |
| 1634 |
enddo |
| 1635 |
enddo |
| 1636 |
enddo |
| 1637 |
C |
| 1638 |
DO 30 I=1,LENC |
| 1639 |
TX1(I) = PRJ(I,K+1) * POI(I,K) |
| 1640 |
QS1(I) = ALF(I,K) + BET(I,K)*POI(I,K) |
| 1641 |
QOL(I,K) = MIN(QS1(I)*RHMAX,QOI(I,K)) |
| 1642 |
|
| 1643 |
HOL(I,K) = TX1(I)*CP + QOL(I,K)*ALHL |
| 1644 |
ETA(I,K) = ZERO |
| 1645 |
TX2(I) = (PRJ(I,K+1) - PRJ(I,K)) * POI(I,K) * CP |
| 1646 |
30 CONTINUE |
| 1647 |
C |
| 1648 |
IF (IC .LT. KM1) THEN |
| 1649 |
DO 3703 L=KM1,IC1,-1 |
| 1650 |
DO 50 I=1,LENC |
| 1651 |
QS1(I) = ALF(I,L) + BET(I,L)*POI(I,L) |
| 1652 |
QOL(I,L) = MIN(QS1(I)*RHMAX,QOI(I,L)) |
| 1653 |
C |
| 1654 |
TEM1 = TX2(I) + PRJ(I,L+1) * POI(I,L) * CP |
| 1655 |
HOL(I,L) = TEM1 + QOL(I,L )* ALHL |
| 1656 |
HST(I,L) = TEM1 + QS1(I) * ALHL |
| 1657 |
|
| 1658 |
TX1(I) = (PRJ(I,L+1) - PRJ(I,L)) * POI(I,L) |
| 1659 |
ETA(I,L) = ETA(I,L+1) + TX1(I)*CPBG |
| 1660 |
TX2(I) = TX2(I) + TX1(I)*CP |
| 1661 |
50 CONTINUE |
| 1662 |
C |
| 1663 |
3703 CONTINUE |
| 1664 |
ENDIF |
| 1665 |
|
| 1666 |
|
| 1667 |
DO 70 I=1,LENC |
| 1668 |
HOL(I,IC) = TX2(I) |
| 1669 |
QS1(I) = ALF(I,IC) + BET(I,IC)*POI(I,IC) |
| 1670 |
QOL(I,IC) = MIN(QS1(I)*RHMAX,QOI(I,IC)) |
| 1671 |
c |
| 1672 |
TEM1 = TX2(I) + PRJ(I,IC1) * POI(I,IC) * CP |
| 1673 |
HOL(I,IC) = TEM1 + QOL(I,IC) * ALHL |
| 1674 |
HST(I,IC) = TEM1 + QS1(I) * ALHL |
| 1675 |
C |
| 1676 |
TX3(I ) = (PRJ(I,IC1) - PRH(I,IC)) * POI(I,IC) |
| 1677 |
ETA(I,IC) = ETA(I,IC1) + CPBG * TX3(I) |
| 1678 |
70 CONTINUE |
| 1679 |
C |
| 1680 |
DO 130 I=1,LENC |
| 1681 |
TX2(I) = HOL(I,K) - HST(I,IC) |
| 1682 |
TX1(I) = ZERO |
| 1683 |
|
| 1684 |
130 CONTINUE |
| 1685 |
C |
| 1686 |
C ENTRAINMENT PARAMETER |
| 1687 |
C |
| 1688 |
DO 160 L=IC,KM1 |
| 1689 |
DO 160 I=1,LENC |
| 1690 |
TX1(I) = TX1(I) + (HST(I,IC) - HOL(I,L)) * (ETA(I,L) - ETA(I,L+1)) |
| 1691 |
160 CONTINUE |
| 1692 |
C |
| 1693 |
LEN1 = 0 |
| 1694 |
LEN2 = 0 |
| 1695 |
ISAV = 0 |
| 1696 |
DO 195 I=1,LENC |
| 1697 |
IF (TX1(I) .GT. ZERO .AND. TX2(I) .GT. ZERO |
| 1698 |
. .AND. rhfrac(i).ne.0.0 ) THEN |
| 1699 |
LEN1 = LEN1 + 1 |
| 1700 |
IA(LEN1) = I |
| 1701 |
ALM(LEN1) = TX2(I) / TX1(I) |
| 1702 |
ENDIF |
| 1703 |
195 CONTINUE |
| 1704 |
C |
| 1705 |
LEN2 = LEN1 |
| 1706 |
if (IC1 .lt. K) then |
| 1707 |
DO 196 I=1,LENC |
| 1708 |
IF (TX2(I) .LE. 0.0 .AND. (HOL(I,K) .GT. HST(I,IC1)) |
| 1709 |
. .AND. rhfrac(i).ne.0.0 ) THEN |
| 1710 |
LEN2 = LEN2 + 1 |
| 1711 |
IA(LEN2) = I |
| 1712 |
ALM(LEN2) = 0.0 |
| 1713 |
ENDIF |
| 1714 |
196 CONTINUE |
| 1715 |
endif |
| 1716 |
C |
| 1717 |
IF (LEN2 .EQ. 0) THEN |
| 1718 |
DO 5010 I=1,LENC*K |
| 1719 |
HST(I,1) = 0.0 |
| 1720 |
QOL(I,1) = 0.0 |
| 1721 |
5010 CONTINUE |
| 1722 |
DO 5020 I=1,LENC |
| 1723 |
PCU(I) = 0.0 |
| 1724 |
5020 CONTINUE |
| 1725 |
RETURN |
| 1726 |
ENDIF |
| 1727 |
LEN11 = LEN1 + 1 |
| 1728 |
C |
| 1729 |
C NORMALIZED MASSFLUX |
| 1730 |
C |
| 1731 |
DO 250 I=1,LEN2 |
| 1732 |
ETA(I,K) = 1.0 |
| 1733 |
II = IA(I) |
| 1734 |
TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1)) |
| 1735 |
TX4(I) = PRS(II,K) |
| 1736 |
250 CONTINUE |
| 1737 |
C |
| 1738 |
DO 252 I=LEN11,LEN2 |
| 1739 |
WFN(I) = 0.0 |
| 1740 |
II = IA(I) |
| 1741 |
IF (HST(II,IC1) .LT. HST(II,IC)) THEN |
| 1742 |
TX6(I) = (HST(II,IC1)-HOL(II,K))/(HST(II,IC1)-HST(II,IC)) |
| 1743 |
ELSE |
| 1744 |
TX6(I) = 0.0 |
| 1745 |
ENDIF |
| 1746 |
TX2(I) = 0.5 * (PRS(II,IC1)+PRS(II,IC1+1)) * (1.0-TX6(I)) |
| 1747 |
* + TX2(I) * TX6(I) |
| 1748 |
252 CONTINUE |
| 1749 |
C |
| 1750 |
CALL ACRITN(LEN2, TX2, TX4, TX3) |
| 1751 |
C |
| 1752 |
DO 260 L=KM1,IC,-1 |
| 1753 |
DO 255 I=1,LEN2 |
| 1754 |
TX1(I) = ETA(IA(I),L) |
| 1755 |
255 CONTINUE |
| 1756 |
DO 260 I=1,LEN2 |
| 1757 |
ETA(I,L) = 1.0 + ALM(I) * TX1(I) |
| 1758 |
260 CONTINUE |
| 1759 |
C |
| 1760 |
C CLOUD WORKFUNCTION |
| 1761 |
C |
| 1762 |
IF (LEN1 .GT. 0) THEN |
| 1763 |
DO 270 I=1,LEN1 |
| 1764 |
II = IA(I) |
| 1765 |
WFN(I) = - GAM(II,IC) * (PRJ(II,IC1) - PRH(II,IC)) |
| 1766 |
* * HST(II,IC) * ETA(I,IC1) |
| 1767 |
270 CONTINUE |
| 1768 |
ENDIF |
| 1769 |
C |
| 1770 |
DO 290 I=1,LEN2 |
| 1771 |
II = IA(I) |
| 1772 |
TX1(I) = HOL(II,K) |
| 1773 |
290 CONTINUE |
| 1774 |
C |
| 1775 |
IF (IC1 .LE. KM1) THEN |
| 1776 |
|
| 1777 |
DO 380 L=KM1,IC1,-1 |
| 1778 |
DO 380 I=1,LEN2 |
| 1779 |
II = IA(I) |
| 1780 |
TEM = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * HOL(II,L) |
| 1781 |
C |
| 1782 |
PCU(I) = PRJ(II,L+1) - PRH(II,L) |
| 1783 |
TEM1 = ETA(I,L+1) * PCU(I) |
| 1784 |
TX1(I) = TX1(I)*PCU(I) |
| 1785 |
C |
| 1786 |
PCU(I) = PRH(II,L) - PRJ(II,L) |
| 1787 |
TEM1 = (TEM1 + ETA(I,L) * PCU(I)) * HST(II,L) |
| 1788 |
TX1(I) = TX1(I) + TEM*PCU(I) |
| 1789 |
C |
| 1790 |
WFN(I) = WFN(I) + (TX1(I) - TEM1) * GAM(II,L) |
| 1791 |
TX1(I) = TEM |
| 1792 |
380 CONTINUE |
| 1793 |
ENDIF |
| 1794 |
C |
| 1795 |
LENA = 0 |
| 1796 |
IF (LEN1 .GT. 0) THEN |
| 1797 |
DO 512 I=1,LEN1 |
| 1798 |
II = IA(I) |
| 1799 |
WFN(I) = WFN(I) + TX1(I) * GAM(II,IC)*(PRJ(II,IC1)-PRH(II,IC)) |
| 1800 |
* - TX3(I) |
| 1801 |
IF (WFN(I) .GT. 0.0) THEN |
| 1802 |
LENA = LENA + 1 |
| 1803 |
I1(LENA) = IA(I) |
| 1804 |
I2(LENA) = I |
| 1805 |
TX1(LENA) = WFN(I) |
| 1806 |
TX2(LENA) = QS1(IA(I)) |
| 1807 |
TX6(LENA) = 1.0 |
| 1808 |
ENDIF |
| 1809 |
512 CONTINUE |
| 1810 |
ENDIF |
| 1811 |
LENB = LENA |
| 1812 |
DO 515 I=LEN11,LEN2 |
| 1813 |
WFN(I) = WFN(I) - TX3(I) |
| 1814 |
IF (WFN(I) .GT. 0.0 .AND. TX6(I) .GT. 0.0) THEN |
| 1815 |
LENB = LENB + 1 |
| 1816 |
I1(LENB) = IA(I) |
| 1817 |
I2(LENB) = I |
| 1818 |
TX1(LENB) = WFN(I) |
| 1819 |
TX2(LENB) = QS1(IA(I)) |
| 1820 |
TX4(LENB) = TX6(I) |
| 1821 |
ENDIF |
| 1822 |
515 CONTINUE |
| 1823 |
C |
| 1824 |
IF (LENB .LE. 0) THEN |
| 1825 |
DO 5030 I=1,LENC*K |
| 1826 |
HST(I,1) = 0.0 |
| 1827 |
QOL(I,1) = 0.0 |
| 1828 |
5030 CONTINUE |
| 1829 |
DO 5040 I=1,LENC |
| 1830 |
PCU(I) = 0.0 |
| 1831 |
5040 CONTINUE |
| 1832 |
RETURN |
| 1833 |
ENDIF |
| 1834 |
|
| 1835 |
C |
| 1836 |
DO 516 I=1,LENB |
| 1837 |
WFN(I) = TX1(I) |
| 1838 |
QS1(I) = TX2(I) |
| 1839 |
516 CONTINUE |
| 1840 |
C |
| 1841 |
DO 520 L=IC,K |
| 1842 |
DO 517 I=1,LENB |
| 1843 |
TX1(I) = ETA(I2(I),L) |
| 1844 |
517 CONTINUE |
| 1845 |
DO 520 I=1,LENB |
| 1846 |
ETA(I,L) = TX1(I) |
| 1847 |
520 CONTINUE |
| 1848 |
C |
| 1849 |
LENA1 = LENA + 1 |
| 1850 |
C |
| 1851 |
DO 510 I=1,LENA |
| 1852 |
II = I1(I) |
| 1853 |
TX8(I) = HST(II,IC) - HOL(II,IC) |
| 1854 |
510 CONTINUE |
| 1855 |
DO 530 I=LENA1,LENB |
| 1856 |
II = I1(I) |
| 1857 |
TX6(I) = TX4(I) |
| 1858 |
TEM = TX6(I) * (HOL(II,IC)-HOL(II,IC1)) + HOL(II,IC1) |
| 1859 |
TX8(I) = HOL(II,K) - TEM |
| 1860 |
|
| 1861 |
TEM1 = TX6(I) * (QOL(II,IC)-QOL(II,IC1)) + QOL(II,IC1) |
| 1862 |
TX5(I) = TEM - TEM1 * ALHL |
| 1863 |
QS1(I) = TEM1 + TX8(I)*(ONE/ALHL) |
| 1864 |
TX3(I) = HOL(II,IC) |
| 1865 |
530 CONTINUE |
| 1866 |
C |
| 1867 |
C |
| 1868 |
DO 620 I=1,LENB |
| 1869 |
II = I1(I) |
| 1870 |
WLQ(I) = QOL(II,K) - QS1(I) * ETA(I,IC) |
| 1871 |
TX7(I) = HOL(II,K) |
| 1872 |
620 CONTINUE |
| 1873 |
DO NT=1,Ntracer |
| 1874 |
DO 621 I=1,LENB |
| 1875 |
II = I1(I) |
| 1876 |
UHT(I,NT) = UOI(II,K+nltop-1,NT)-UOI(II,IC+nltop-1,NT) * ETA(I,IC) |
| 1877 |
621 CONTINUE |
| 1878 |
ENDDO |
| 1879 |
C |
| 1880 |
DO 635 L=KM1,IC,-1 |
| 1881 |
DO 630 I=1,LENB |
| 1882 |
II = I1(I) |
| 1883 |
TEM = ETA(I,L) - ETA(I,L+1) |
| 1884 |
WLQ(I) = WLQ(I) + TEM * QOL(II,L) |
| 1885 |
630 CONTINUE |
| 1886 |
635 CONTINUE |
| 1887 |
DO NT=1,Ntracer |
| 1888 |
DO L=KM1,IC,-1 |
| 1889 |
DO I=1,LENB |
| 1890 |
II = I1(I) |
| 1891 |
TEM = ETA(I,L) - ETA(I,L+1) |
| 1892 |
UHT(I,NT) = UHT(I,NT) + TEM * UOI(II,L+nltop-1,NT) |
| 1893 |
ENDDO |
| 1894 |
ENDDO |
| 1895 |
ENDDO |
| 1896 |
C |
| 1897 |
C CALCULATE GS AND PART OF AKM (THAT REQUIRES ETA) |
| 1898 |
C |
| 1899 |
DO 690 I=1,LENB |
| 1900 |
II = I1(I) |
| 1901 |
c TX7(I) = HOL(II,K) |
| 1902 |
TEM = (POI(II,KM1) - POI(II,K)) / (PRH(II,K) - PRH(II,KM1)) |
| 1903 |
HOL(I,K) = TEM * (PRJ(II,K)-PRH(II,KM1))*PRH(II,K)*PRI(II,K) |
| 1904 |
HOL(I,KM1) = TEM * (PRH(II,K)-PRJ(II,K))*PRH(II,KM1)*PRI(II,KM1) |
| 1905 |
AKM(I) = ZERO |
| 1906 |
TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1)) |
| 1907 |
690 CONTINUE |
| 1908 |
|
| 1909 |
IF (IC1 .LE. KM1) THEN |
| 1910 |
DO 750 L=KM1,IC1,-1 |
| 1911 |
DO 750 I=1,LENB |
| 1912 |
II = I1(I) |
| 1913 |
TEM = (POI(II,L-1) - POI(II,L)) * ETA(I,L) |
| 1914 |
* / (PRH(II,L) - PRH(II,L-1)) |
| 1915 |
C |
| 1916 |
HOL(I,L) = TEM * (PRJ(II,L)-PRH(II,L-1)) * PRH(II,L) |
| 1917 |
* * PRI(II,L) + HOL(I,L) |
| 1918 |
HOL(I,L-1) = TEM * (PRH(II,L)-PRJ(II,L)) * PRH(II,L-1) |
| 1919 |
* * PRI(II,L-1) |
| 1920 |
C |
| 1921 |
AKM(I) = AKM(I) - HOL(I,L) |
| 1922 |
* * (ETA(I,L) * (PRH(II,L)-PRJ(II,L)) + |
| 1923 |
* ETA(I,L+1) * (PRJ(II,L+1)-PRH(II,L))) / PRH(II,L) |
| 1924 |
750 CONTINUE |
| 1925 |
ENDIF |
| 1926 |
C |
| 1927 |
C |
| 1928 |
CALL RNCL(LENB, TX2, TX1, CLF) |
| 1929 |
|
| 1930 |
DO 770 I=1,LENB |
| 1931 |
TX2(I) = (ONE - TX1(I)) * WLQ(I) |
| 1932 |
WLQ(I) = TX1(I) * WLQ(I) |
| 1933 |
C |
| 1934 |
TX1(I) = HOL(I,IC) |
| 1935 |
770 CONTINUE |
| 1936 |
DO 790 I=LENA1, LENB |
| 1937 |
II = I1(I) |
| 1938 |
TX1(I) = TX1(I) + (TX5(I)-TX3(I)+QOL(II,IC)*ALHL)*(PRI(II,IC)/CP) |
| 1939 |
790 CONTINUE |
| 1940 |
|
| 1941 |
DO 800 I=1,LENB |
| 1942 |
HOL(I,IC) = TX1(I) - TX2(I) * ALBCP * PRI(I1(I),IC) |
| 1943 |
800 CONTINUE |
| 1944 |
|
| 1945 |
IF (LENA .GT. 0) THEN |
| 1946 |
DO 810 I=1,LENA |
| 1947 |
II = I1(I) |
| 1948 |
AKM(I) = AKM(I) - ETA(I,IC1) * (PRJ(II,IC1) - PRH(II,IC)) |
| 1949 |
* * TX1(I) / PRH(II,IC) |
| 1950 |
810 CONTINUE |
| 1951 |
ENDIF |
| 1952 |
c |
| 1953 |
C CALCULATE GH |
| 1954 |
C |
| 1955 |
DO 830 I=1,LENB |
| 1956 |
II = I1(I) |
| 1957 |
TX3(I) = QOL(II,KM1) - QOL(II,K) |
| 1958 |
GMH(I,K) = HOL(I,K) + TX3(I) * PRI(II,K) * (ALBCP) |
| 1959 |
|
| 1960 |
AKM(I) = AKM(I) + GAM(II,KM1)*(PRJ(II,K)-PRH(II,KM1)) |
| 1961 |
* * GMH(I,K) |
| 1962 |
TX3(I) = zero |
| 1963 |
830 CONTINUE |
| 1964 |
C |
| 1965 |
IF (IC1 .LE. KM1) THEN |
| 1966 |
DO 840 L=KM1,IC1,-1 |
| 1967 |
DO 840 I=1,LENB |
| 1968 |
II = I1(I) |
| 1969 |
TX2(I) = TX3(I) |
| 1970 |
TX3(I) = (QOL(II,L-1) - QOL(II,L)) * ETA(I,L) |
| 1971 |
TX2(I) = TX2(I) + TX3(I) |
| 1972 |
C |
| 1973 |
GMH(I,L) = HOL(I,L) + TX2(I) * PRI(II,L) * (ALBCP*HALF) |
| 1974 |
840 CONTINUE |
| 1975 |
C |
| 1976 |
C |
| 1977 |
ENDIF |
| 1978 |
DO 850 I=LENA1,LENB |
| 1979 |
TX3(I) = TX3(I) + TWOBAL |
| 1980 |
* * (TX7(I) - TX8(I) - TX5(I) - QOL(I1(I),IC)*ALHL) |
| 1981 |
850 CONTINUE |
| 1982 |
DO 860 I=1,LENB |
| 1983 |
GMH(I,IC) = TX1(I) + PRI(I1(I),IC) * ONEBCP |
| 1984 |
* * (TX3(I)*(ALHL*HALF) + ETA(I,IC) * TX8(I)) |
| 1985 |
860 CONTINUE |
| 1986 |
C |
| 1987 |
C CALCULATE HC PART OF AKM |
| 1988 |
C |
| 1989 |
IF (IC1 .LE. KM1) THEN |
| 1990 |
DO 870 I=1,LENB |
| 1991 |
TX1(I) = GMH(I,K) |
| 1992 |
870 CONTINUE |
| 1993 |
DO 3725 L=KM1,IC1,-1 |
| 1994 |
DO 880 I=1,LENB |
| 1995 |
II = I1(I) |
| 1996 |
TX1(I) = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * GMH(I,L) |
| 1997 |
TX2(I) = GAM(II,L-1) * (PRJ(II,L) - PRH(II,L-1)) |
| 1998 |
880 CONTINUE |
| 1999 |
C |
| 2000 |
IF (L .EQ. IC1) THEN |
| 2001 |
DO 890 I=LENA1,LENB |
| 2002 |
TX2(I) = ZERO |
| 2003 |
890 CONTINUE |
| 2004 |
ENDIF |
| 2005 |
DO 900 I=1,LENB |
| 2006 |
II = I1(I) |
| 2007 |
AKM(I) = AKM(I) + TX1(I) * |
| 2008 |
* (TX2(I) + GAM(II,L)*(PRH(II,L)-PRJ(II,L))) |
| 2009 |
900 CONTINUE |
| 2010 |
3725 CONTINUE |
| 2011 |
ENDIF |
| 2012 |
C |
| 2013 |
DO 920 I=LENA1,LENB |
| 2014 |
II = I1(I) |
| 2015 |
TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1)) |
| 2016 |
* + 0.5*(PRS(II,IC+2) - PRS(II,IC)) * (ONE-TX6(I)) |
| 2017 |
c |
| 2018 |
TX1(I) = PRS(II,IC1) |
| 2019 |
TX5(I) = 0.5 * (PRS(II,IC1) + PRS(II,IC+2)) |
| 2020 |
C |
| 2021 |
IF ((TX2(I) .GE. TX1(I)) .AND. (TX2(I) .LT. TX5(I))) THEN |
| 2022 |
TX6(I) = ONE - (TX2(I) - TX1(I)) / (TX5(I) - TX1(I)) |
| 2023 |
C |
| 2024 |
TEM = PRI(II,IC1) / PRI(II,IC) |
| 2025 |
HOL(I,IC1) = HOL(I,IC1) + HOL(I,IC) * TEM |
| 2026 |
HOL(I,IC) = ZERO |
| 2027 |
C |
| 2028 |
GMH(I,IC1) = GMH(I,IC1) + GMH(I,IC) * TEM |
| 2029 |
GMH(I,IC) = ZERO |
| 2030 |
ELSEIF (TX2(I) .LT. TX1(I)) THEN |
| 2031 |
TX6(I) = 1.0 |
| 2032 |
ELSE |
| 2033 |
TX6(I) = 0.0 |
| 2034 |
ENDIF |
| 2035 |
920 CONTINUE |
| 2036 |
C |
| 2037 |
C |
| 2038 |
DO I=1,LENC |
| 2039 |
PCU(I) = 0.0 |
| 2040 |
ENDDO |
| 2041 |
|
| 2042 |
DO 970 I=1,LENB |
| 2043 |
II = I1(I) |
| 2044 |
IF (AKM(I) .LT. ZERO .AND. WLQ(I) .GE. 0.0) THEN |
| 2045 |
WFN(I) = - TX6(I) * WFN(I) * RASALF / AKM(I) |
| 2046 |
ELSE |
| 2047 |
WFN(I) = ZERO |
| 2048 |
ENDIF |
| 2049 |
TEM = (PRS(II,K+1)-PRS(II,K))*(CMB2PA*FRAC) |
| 2050 |
WFN(I) = MIN(WFN(I), TEM) |
| 2051 |
C |
| 2052 |
C compute cloud amount |
| 2053 |
C |
| 2054 |
CC TX1(I) = CLN(II) |
| 2055 |
CC IF (WFN(I) .GT. CRTMSF) TX1(I) = TX1(I) + CLF(I) |
| 2056 |
CC IF (TX1(I) .GT. ONE) TX1(I) = ONE |
| 2057 |
C |
| 2058 |
C PRECIPITATION |
| 2059 |
C |
| 2060 |
PCU(II) = WLQ(I) * WFN(I) * ONEBG |
| 2061 |
C |
| 2062 |
C CUMULUS FRICTION AT THE BOTTOM LAYER |
| 2063 |
C |
| 2064 |
TX4(I) = WFN(I) * (1.0/ALHL) |
| 2065 |
TX5(I) = WFN(I) * ONEBCP |
| 2066 |
970 CONTINUE |
| 2067 |
C |
| 2068 |
C compute cloud mass flux for diagnostic output |
| 2069 |
C |
| 2070 |
DO L = IC,K |
| 2071 |
DO I=1,LENB |
| 2072 |
II = I1(I) |
| 2073 |
if(L.lt.K)then |
| 2074 |
CMASS(II,L) = ETA(I,L+1) * WFN(I) * ONEBG |
| 2075 |
else |
| 2076 |
CMASS(II,L) = WFN(I) * ONEBG |
| 2077 |
endif |
| 2078 |
ENDDO |
| 2079 |
ENDDO |
| 2080 |
C |
| 2081 |
CC DO 975 I=1,LENB |
| 2082 |
CC II = I1(I) |
| 2083 |
CC CLN(II) = TX1(I) |
| 2084 |
CC975 CONTINUE |
| 2085 |
C |
| 2086 |
C THETA AND Q CHANGE DUE TO CLOUD TYPE IC |
| 2087 |
C |
| 2088 |
|
| 2089 |
c TEMA = 0.0 |
| 2090 |
c TEMB = 0.0 |
| 2091 |
DO 990 L=IC,K |
| 2092 |
DO 980 I=1,LENB |
| 2093 |
II = I1(I) |
| 2094 |
TEM = (GMH(I,L) - HOL(I,L)) * TX4(I) |
| 2095 |
TEM1 = HOL(I,L) * TX5(I) |
| 2096 |
C |
| 2097 |
TCU(II,L) = TEM1 / PRH(II,L) |
| 2098 |
QCU(II,L) = TEM |
| 2099 |
980 CONTINUE |
| 2100 |
|
| 2101 |
c I = I1(IP1) |
| 2102 |
c |
| 2103 |
c TEM = (PRS(I,L+1)-PRS(I,L)) * (ONEBG*100.0) |
| 2104 |
c TEMA = TEMA + TCU(I,L) * PRH(I,L) * TEM * (CP/ALHL) |
| 2105 |
c TEMB = TEMB + QCU(I,L) * TEM |
| 2106 |
C |
| 2107 |
990 CONTINUE |
| 2108 |
C |
| 2109 |
c Compute Tracer Tendencies |
| 2110 |
c ------------------------- |
| 2111 |
do nt = 1,ntracer |
| 2112 |
|
| 2113 |
c Tracer Tendency at the Bottom Layer |
| 2114 |
c ----------------------------------- |
| 2115 |
DO 995 I=1,LENB |
| 2116 |
II = I1(I) |
| 2117 |
TEM = half*TX5(I) * PRI(II,K) |
| 2118 |
TX1(I) = (UOI(II,KM1+nltop-1,nt) - UOI(II,K+nltop-1,nt)) |
| 2119 |
ucu(II,K,nt) = TEM * TX1(I) |
| 2120 |
995 CONTINUE |
| 2121 |
|
| 2122 |
c Tracer Tendency at all other Levels |
| 2123 |
c ----------------------------------- |
| 2124 |
DO 1020 L=KM1,IC1,-1 |
| 2125 |
DO 1010 I=1,LENB |
| 2126 |
II = I1(I) |
| 2127 |
TEM = half*TX5(I) * PRI(II,L) |
| 2128 |
TEM1 = TX1(I) |
| 2129 |
TX1(I) = (UOI(II,L-1+nltop-1,nt)-UOI(II,L+nltop-1,nt)) * ETA(I,L) |
| 2130 |
TX3(I) = (TX1(I) + TEM1) * TEM |
| 2131 |
1010 CONTINUE |
| 2132 |
DO 1020 I=1,LENB |
| 2133 |
II = I1(I) |
| 2134 |
ucu(II,L,nt) = TX3(I) |
| 2135 |
1020 CONTINUE |
| 2136 |
|
| 2137 |
DO 1030 I=1,LENB |
| 2138 |
II = I1(I) |
| 2139 |
IF (TX6(I) .GE. 1.0) THEN |
| 2140 |
TEM = half*TX5(I) * PRI(II,IC) |
| 2141 |
ELSE |
| 2142 |
TEM = 0.0 |
| 2143 |
ENDIF |
| 2144 |
TX1(I) = (TX1(I) + UHT(I,nt) + UHT(I,nt)) * TEM |
| 2145 |
1030 CONTINUE |
| 2146 |
DO 1040 I=1,LENB |
| 2147 |
II = I1(I) |
| 2148 |
ucu(II,IC,nt) = TX1(I) |
| 2149 |
1040 CONTINUE |
| 2150 |
|
| 2151 |
enddo |
| 2152 |
C |
| 2153 |
C PENETRATIVE CONVECTION CALCULATION OVER |
| 2154 |
C |
| 2155 |
|
| 2156 |
RETURN |
| 2157 |
END |
| 2158 |
SUBROUTINE RNCL(LEN, PL, RNO, CLF) |
| 2159 |
C |
| 2160 |
C |
| 2161 |
C********************************************************************* |
| 2162 |
C********************** Relaxed Arakawa-Schubert ********************* |
| 2163 |
C************************ SUBROUTINE RNCL ************************ |
| 2164 |
C**************************** 23 July 1992 *************************** |
| 2165 |
C********************************************************************* |
| 2166 |
|
| 2167 |
PARAMETER (P5=500.0, P8=800.0, PT8=0.8, PT2=0.2) |
| 2168 |
PARAMETER (PFAC=PT2/(P8-P5)) |
| 2169 |
C |
| 2170 |
PARAMETER (P4=400.0, P6=401.0) |
| 2171 |
PARAMETER (P7=700.0, P9=900.0) |
| 2172 |
PARAMETER (CUCLD=0.5,CFAC=CUCLD/(P6-P4)) |
| 2173 |
C |
| 2174 |
DIMENSION PL(LEN), RNO(LEN), CLF(LEN) |
| 2175 |
|
| 2176 |
DO 10 I=1,LEN |
| 2177 |
rno(i) = 1.0 |
| 2178 |
ccc if( pl(i).le.400.0 ) rno(i) = max( 0.75, 1.0-0.0025*(400.0-pl(i)) ) |
| 2179 |
|
| 2180 |
ccc IF ( PL(I).GE.P7 .AND. PL(I).LE.P9 ) THEN |
| 2181 |
ccc RNO(I) = ((P9-PL(I))/(P9-P7)) **2 |
| 2182 |
ccc ELSE IF (PL(I).GT.P9) THEN |
| 2183 |
ccc RNO(I) = 0. |
| 2184 |
ccc ENDIF |
| 2185 |
|
| 2186 |
CLF(I) = CUCLD |
| 2187 |
C |
| 2188 |
CARIESIF (PL(I) .GE. P5 .AND. PL(I) .LE. P8) THEN |
| 2189 |
CARIES RNO(I) = (P8-PL(I))*PFAC + PT8 |
| 2190 |
CARIESELSEIF (PL(I) .GT. P8 ) THEN |
| 2191 |
CARIES RNO(I) = PT8 |
| 2192 |
CARIESENDIF |
| 2193 |
CARIES |
| 2194 |
IF (PL(I) .GE. P4 .AND. PL(I) .LE. P6) THEN |
| 2195 |
CLF(I) = (P6-PL(I))*CFAC |
| 2196 |
ELSEIF (PL(I) .GT. P6 ) THEN |
| 2197 |
CLF(I) = 0.0 |
| 2198 |
ENDIF |
| 2199 |
10 CONTINUE |
| 2200 |
C |
| 2201 |
RETURN |
| 2202 |
END |
| 2203 |
SUBROUTINE ACRITN ( LEN,PL,PLB,ACR ) |
| 2204 |
|
| 2205 |
C********************************************************************* |
| 2206 |
C********************** Relaxed Arakawa-Schubert ********************* |
| 2207 |
C************************** SUBROUTINE ACRIT ********************* |
| 2208 |
C****************** modified August 28, 1996 L.Takacs ************ |
| 2209 |
C**** ***** |
| 2210 |
C**** Note: Data obtained from January Mean After-Analysis ***** |
| 2211 |
C**** from 4x5 46-layer GEOS Assimilation ***** |
| 2212 |
C**** ***** |
| 2213 |
C********************************************************************* |
| 2214 |
|
| 2215 |
real PL(LEN), PLB(LEN), ACR(LEN) |
| 2216 |
|
| 2217 |
parameter (lma=18) |
| 2218 |
real p(lma) |
| 2219 |
real a(lma) |
| 2220 |
|
| 2221 |
data p / 93.81, 111.65, 133.46, 157.80, 186.51, |
| 2222 |
. 219.88, 257.40, 301.21, 352.49, 409.76, |
| 2223 |
. 471.59, 535.04, 603.33, 672.79, 741.12, |
| 2224 |
. 812.52, 875.31, 930.20/ |
| 2225 |
|
| 2226 |
data a / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364, |
| 2227 |
. 1.01971, .65846, .45867, .38687, .31002, |
| 2228 |
. .25574, .20347, .17254, .15260, .16756, |
| 2229 |
. .09916, .10360, .05880/ |
| 2230 |
|
| 2231 |
|
| 2232 |
do L=1,lma-1 |
| 2233 |
do i=1,len |
| 2234 |
if( pl(i).ge.p(L) .and. |
| 2235 |
. pl(i).le.p(L+1)) then |
| 2236 |
temp = ( pl(i)-p(L) )/( p(L+1)-p(L) ) |
| 2237 |
acr(i) = a(L+1)*temp + a(L)*(1-temp) |
| 2238 |
endif |
| 2239 |
enddo |
| 2240 |
enddo |
| 2241 |
|
| 2242 |
do i=1,len |
| 2243 |
if( pl(i).lt.p(1) ) acr(i) = a(1) |
| 2244 |
if( pl(i).gt.p(lma) ) acr(i) = a(lma) |
| 2245 |
enddo |
| 2246 |
|
| 2247 |
do i=1,len |
| 2248 |
acr(i) = acr(i) * (plb(i)-pl(i)) |
| 2249 |
enddo |
| 2250 |
|
| 2251 |
RETURN |
| 2252 |
END |
| 2253 |
SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DSIG,PLKE, |
| 2254 |
1 PLK,TH,TEMP1,TEMP2,TEMP3,ITMP1,ITMP2,RCON,RLAR,CLSBTH,tmscl, |
| 2255 |
2 tmfrc,cp,gravity,alhl,gamfac,cldlz,RHCRIT,offset,alpha) |
| 2256 |
|
| 2257 |
PARAMETER (ZM1P04 = -1.04E-4 ) |
| 2258 |
PARAMETER (ZERO = 0.) |
| 2259 |
PARAMETER (TWO89= 2.89E-5) |
| 2260 |
PARAMETER ( ZP44= 0.44) |
| 2261 |
PARAMETER ( ZP01= 0.01) |
| 2262 |
PARAMETER ( ZP1 = 0.1 ) |
| 2263 |
PARAMETER ( ZP001= 0.001) |
| 2264 |
PARAMETER ( HALF= 0.5) |
| 2265 |
PARAMETER ( ZP578 = 0.578 ) |
| 2266 |
PARAMETER ( ONE = 1.) |
| 2267 |
PARAMETER ( THOUSAND = 1000.) |
| 2268 |
PARAMETER ( Z3600 = 3600.) |
| 2269 |
C |
| 2270 |
DIMENSION TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY), |
| 2271 |
$ PL(IRUN,NLAY),CLFRAC(IRUN,NLAY),SP(IRUN),TEMP1(IRUN,NLAY), |
| 2272 |
$ TEMP2(IRUN,NLAY),PLKE(IRUN,NLAY), |
| 2273 |
$ RCON(IRUN),RLAR(IRUN),DSIG(NLAY),PLK(IRUN,NLAY),TH(IRUN,NLAY), |
| 2274 |
$ TEMP3(IRUN,NLAY),ITMP1(IRUN,NLAY), |
| 2275 |
$ ITMP2(IRUN,NLAY),CLSBTH(IRUN,NLAY) |
| 2276 |
C |
| 2277 |
DIMENSION EVP9(IRUN,NLAY) |
| 2278 |
real water(irun),crystal(irun) |
| 2279 |
real watevap(irun),iceevap(irun) |
| 2280 |
real fracwat,fracice, tice,rh,fact,dum |
| 2281 |
|
| 2282 |
real cldlz(irun,nlay) |
| 2283 |
real rhcrit(irun,nlay), rainmax(irun) |
| 2284 |
real offset, alpha |
| 2285 |
|
| 2286 |
c Explicit Inline Directives |
| 2287 |
c -------------------------- |
| 2288 |
#if CRAY |
| 2289 |
#if f77 |
| 2290 |
cfpp$ expand (qsat) |
| 2291 |
#endif |
| 2292 |
#endif |
| 2293 |
|
| 2294 |
tice = getcon('FREEZING-POINT') |
| 2295 |
|
| 2296 |
fracwat = 0.70 |
| 2297 |
fracice = 0.01 |
| 2298 |
|
| 2299 |
NLAYM1 = NLAY - 1 |
| 2300 |
IRNLAY = IRUN*NLAY |
| 2301 |
IRNLM1 = IRUN*(NLAY-1) |
| 2302 |
|
| 2303 |
RPHF = Z3600/tmscl |
| 2304 |
|
| 2305 |
ELOCP = alhl/cp |
| 2306 |
CPOG = cp/gravity |
| 2307 |
|
| 2308 |
DO I = 1,IRUN |
| 2309 |
RLAR(I) = 0. |
| 2310 |
water(i) = 0. |
| 2311 |
crystal(i) = 0. |
| 2312 |
ENDDO |
| 2313 |
|
| 2314 |
do L = 1,nlay |
| 2315 |
do i = 1,irun |
| 2316 |
EVP9(i,L) = 0. |
| 2317 |
TEMP1(i,L) = 0. |
| 2318 |
TEMP2(i,L) = 0. |
| 2319 |
TEMP3(i,L) = 0. |
| 2320 |
CLSBTH(i,L) = 0. |
| 2321 |
cldlz(i,L) = 0. |
| 2322 |
enddo |
| 2323 |
enddo |
| 2324 |
|
| 2325 |
C RHO(ZERO) / RHO FOR TERMINAL VELOCITY APPROX. |
| 2326 |
c --------------------------------------------- |
| 2327 |
DO L = 1,NLAY |
| 2328 |
DO I = 1,IRUN |
| 2329 |
TEMP2(I,L) = PL(I,L)*ZP001 |
| 2330 |
TEMP2(I,L) = SQRT(TEMP2(I,L)) |
| 2331 |
ENDDO |
| 2332 |
ENDDO |
| 2333 |
|
| 2334 |
C INVERSE OF MASS IN EACH LAYER |
| 2335 |
c ----------------------------- |
| 2336 |
DO L = 1,NLAY |
| 2337 |
DO I = 1,IRUN |
| 2338 |
TEMP3(I,L) = SP(I) * DSIG(L) |
| 2339 |
TEMP3(I,L) = GRAVITY*ZP01 / TEMP3(I,L) |
| 2340 |
ENDDO |
| 2341 |
ENDDO |
| 2342 |
|
| 2343 |
C DO LOOP FOR MOISTURE EVAPORATION ABILITY AND CONVEC EVAPORATION. |
| 2344 |
c ---------------------------------------------------------------- |
| 2345 |
DO 100 L=1,NLAY |
| 2346 |
|
| 2347 |
DO I = 1,IRUN |
| 2348 |
TEMP1(I,3) = TL(I,L) |
| 2349 |
TEMP1(I,4) = QL(I,L) |
| 2350 |
ENDDO |
| 2351 |
|
| 2352 |
DO 50 N=1,2 |
| 2353 |
IF(N.EQ.1)RELAX=HALF |
| 2354 |
IF(N.GT.1)RELAX=ONE |
| 2355 |
|
| 2356 |
DO I = 1,IRUN |
| 2357 |
call qsat ( temp1(i,3),pl(i,L),temp1(i,2),temp1(i,6),.true. ) |
| 2358 |
TEMP1(I,5)=TEMP1(I,2)-TEMP1(I,4) |
| 2359 |
TEMP1(I,6)=TEMP1(I,6)*ELOCP |
| 2360 |
TEMP1(I,5)=TEMP1(I,5)/(ONE+TEMP1(I,6)) |
| 2361 |
TEMP1(I,4)=TEMP1(I,4)+TEMP1(I,5)*RELAX |
| 2362 |
TEMP1(I,3)=TEMP1(I,3)-TEMP1(I,5)*ELOCP*RELAX |
| 2363 |
ENDDO |
| 2364 |
50 CONTINUE |
| 2365 |
|
| 2366 |
DO I = 1,IRUN |
| 2367 |
EVP9(I,L) = (TEMP1(I,4) - QL(I,L))/TEMP3(I,L) |
| 2368 |
C convective detrained water |
| 2369 |
cldlz(i,L) = rain(i,L)*temp3(i,L) |
| 2370 |
if( tl(i,L).gt.tice-20.) then |
| 2371 |
water(i) = water(i) + rain(i,L) |
| 2372 |
else |
| 2373 |
crystal(i) = crystal(i) + rain(i,L) |
| 2374 |
endif |
| 2375 |
ENDDO |
| 2376 |
|
| 2377 |
C********************************************************************** |
| 2378 |
C FOR CONVECTIVE PRECIP, FIND THE "EVAPORATION EFFICIENCY" USING * |
| 2379 |
C KESSLERS PARAMETERIZATION * |
| 2380 |
C********************************************************************** |
| 2381 |
|
| 2382 |
DO 20 I=1,IRUN |
| 2383 |
|
| 2384 |
iceevap(i) = 0. |
| 2385 |
watevap(i) = 0. |
| 2386 |
|
| 2387 |
if( (evp9(i,L).gt.0.) .and. (crystal(i).gt.0.) ) then |
| 2388 |
iceevap(I) = EVP9(I,L)*fracice |
| 2389 |
IF(iceevap(i).GE.crystal(i)) iceevap(i) = crystal(i) |
| 2390 |
EVP9(I,L)=EVP9(I,L)-iceevap(I) |
| 2391 |
crystal(i) = crystal(i) - iceevap(i) |
| 2392 |
endif |
| 2393 |
|
| 2394 |
C and now warm precipitate |
| 2395 |
if( (evp9(i,L).gt.0.) .and. (water(i).gt.0.) ) then |
| 2396 |
exparg = ZM1P04*tmscl*((water(i)*RPHF*TEMP2(I,L))**ZP578) |
| 2397 |
AREARAT = ONE-(EXP(EXPARG)) |
| 2398 |
watevap(I) = EVP9(I,L)*AREARAT*fracwat |
| 2399 |
IF(watevap(I).GE.water(i)) watevap(I) = water(i) |
| 2400 |
EVP9(I,L)=EVP9(I,L)-watevap(I) |
| 2401 |
water(i) = water(i) - watevap(i) |
| 2402 |
endif |
| 2403 |
|
| 2404 |
QL(I,L) = QL(I,L)+(iceevap(i)+watevap(i))*TEMP3(I,L) |
| 2405 |
TL(I,L) = TL(I,L)-(iceevap(i)+watevap(i))*TEMP3(I,L)*ELOCP |
| 2406 |
|
| 2407 |
20 CONTINUE |
| 2408 |
|
| 2409 |
100 CONTINUE |
| 2410 |
|
| 2411 |
do i = 1,irun |
| 2412 |
rcon(i) = water(i) + crystal(i) |
| 2413 |
enddo |
| 2414 |
|
| 2415 |
C********************************************************************** |
| 2416 |
C Large Scale Precip |
| 2417 |
C********************************************************************** |
| 2418 |
|
| 2419 |
DO 200 L=1,NLAY |
| 2420 |
DO I = 1,IRUN |
| 2421 |
rainmax(i) = rhcrit(i,L)*evp9(i,L) + |
| 2422 |
. ql(i,L)*(rhcrit(i,L)-1.)/temp3(i,L) |
| 2423 |
|
| 2424 |
if (rainmax(i).LE.0.0) then |
| 2425 |
call qsat( tl(i,L),pl(i,L),rh,dum,.false.) |
| 2426 |
rh = ql(i,L)/rh |
| 2427 |
|
| 2428 |
if( rhcrit(i,L).eq.1.0 ) then |
| 2429 |
fact = 1.0 |
| 2430 |
else |
| 2431 |
fact = min( 1.0, alpha + (1.0-alpha)*( rh-rhcrit(i,L)) / |
| 2432 |
1 (1.0-rhcrit(i,L)) ) |
| 2433 |
endif |
| 2434 |
|
| 2435 |
C Do not allow clouds above 10 mb |
| 2436 |
if( pl(i,L).ge.10.0 ) CLSBTH(I,L) = fact |
| 2437 |
RLAR(I) = RLAR(I)-rainmax(I) |
| 2438 |
QL(I,L) = QL(I,L)+rainmax(I)*TEMP3(I,L) |
| 2439 |
TL(I,L) = TL(I,L)-rainmax(I)*TEMP3(I,L)*ELOCP |
| 2440 |
C Large-scale water |
| 2441 |
cldlz(i,L) = cldlz(i,L) - rainmax(i)*temp3(i,L) |
| 2442 |
ENDIF |
| 2443 |
ENDDO |
| 2444 |
|
| 2445 |
DO I=1,IRUN |
| 2446 |
IF((RLAR(I).GT.0.0).AND.(rainmax(I).GT.0.0))THEN |
| 2447 |
RPOW=(RLAR(I)*RPHF*TEMP2(I,L))**ZP578 |
| 2448 |
EXPARG = ZM1P04*tmscl*RPOW |
| 2449 |
AREARAT = ONE-(EXP(EXPARG)) |
| 2450 |
TEMP1(I,7) = rainmax(I)*AREARAT |
| 2451 |
IF(TEMP1(I,7).GE.RLAR(I)) TEMP1(I,7) = RLAR(I) |
| 2452 |
RLAR(I) = RLAR(I)-TEMP1(I,7) |
| 2453 |
QL(I,L) = QL(I,L)+TEMP1(I,7)*TEMP3(I,L) |
| 2454 |
TL(I,L) = TL(I,L)-TEMP1(I,7)*TEMP3(I,L)*ELOCP |
| 2455 |
ENDIF |
| 2456 |
ENDDO |
| 2457 |
|
| 2458 |
200 CONTINUE |
| 2459 |
|
| 2460 |
RETURN |
| 2461 |
END |
| 2462 |
subroutine srclouds (th,q,plk,pl,plke,cloud,cldwat,irun,irise, |
| 2463 |
1 rhc,offset,alpha) |
| 2464 |
C*********************************************************************** |
| 2465 |
C |
| 2466 |
C PURPOSE: |
| 2467 |
C ======== |
| 2468 |
C Compute non-precipitating cloud fractions |
| 2469 |
C based on Slingo and Ritter (1985). |
| 2470 |
C Remove cloudiness where conditionally unstable. |
| 2471 |
C |
| 2472 |
C INPUT: |
| 2473 |
C ====== |
| 2474 |
C th ......... Potential Temperature (irun,irise) |
| 2475 |
C q .......... Specific Humidity (irun,irise) |
| 2476 |
C plk ........ P**Kappa at mid-layer (irun,irise) |
| 2477 |
C pl ......... Pressure at mid-layer (irun,irise) |
| 2478 |
C plke ....... P**Kappa at edge (irun,irise+1) |
| 2479 |
C irun ....... Horizontal dimension |
| 2480 |
C irise ...... Vertical dimension |
| 2481 |
C |
| 2482 |
C OUTPUT: |
| 2483 |
C ======= |
| 2484 |
C cloud ...... Cloud Fraction (irun,irise) |
| 2485 |
C |
| 2486 |
C*********************************************************************** |
| 2487 |
C* GODDARD LABORATORY FOR ATMOSPHERES * |
| 2488 |
C*********************************************************************** |
| 2489 |
|
| 2490 |
implicit none |
| 2491 |
integer irun,irise |
| 2492 |
|
| 2493 |
real th(irun,irise) |
| 2494 |
real q(irun,irise) |
| 2495 |
real plk(irun,irise) |
| 2496 |
real pl(irun,irise) |
| 2497 |
real plke(irun,irise+1) |
| 2498 |
|
| 2499 |
real tempth(irun) |
| 2500 |
real tempqs(irun) |
| 2501 |
real dhstar(irun) |
| 2502 |
real cloud(irun,irise) |
| 2503 |
real cldwat(irun,irise) |
| 2504 |
real qs(irun,irise) |
| 2505 |
|
| 2506 |
real cp, alhl, getcon, akap, pcheck |
| 2507 |
real ratio, temp, pke, elocp |
| 2508 |
real rhcrit,rh,dum,pbar,tbar |
| 2509 |
integer i,L,ntradesu,ntradesl |
| 2510 |
|
| 2511 |
real factor |
| 2512 |
real rhc(irun,irise) |
| 2513 |
real offset,alpha |
| 2514 |
|
| 2515 |
c Explicit Inline Directives |
| 2516 |
c -------------------------- |
| 2517 |
#if CRAY |
| 2518 |
#if f77 |
| 2519 |
cfpp$ expand (qsat) |
| 2520 |
#endif |
| 2521 |
#endif |
| 2522 |
|
| 2523 |
cp = getcon('CP') |
| 2524 |
alhl = getcon('LATENT HEAT COND') |
| 2525 |
elocp = alhl/cp |
| 2526 |
akap = getcon('KAPPA') |
| 2527 |
|
| 2528 |
do L = 1,irise |
| 2529 |
do i = 1,irun |
| 2530 |
temp = th(i,L)*plk(i,L) |
| 2531 |
call qsat ( temp,pl(i,L),qs(i,L),dum,.false. ) |
| 2532 |
enddo |
| 2533 |
enddo |
| 2534 |
|
| 2535 |
do L = 2,irise |
| 2536 |
do i = 1,irun |
| 2537 |
rh = q(i,L)/qs(i,L) |
| 2538 |
|
| 2539 |
rhcrit = rhc(i,L) - offset |
| 2540 |
ratio = alpha*(rh-rhcrit)/offset |
| 2541 |
|
| 2542 |
if(cloud(i,L).eq. 0.0 .and. ratio.gt.0.0 ) then |
| 2543 |
cloud(i,L) = min( ratio,1.0 ) |
| 2544 |
endif |
| 2545 |
|
| 2546 |
enddo |
| 2547 |
enddo |
| 2548 |
|
| 2549 |
c Reduce clouds from conditionally unstable layer |
| 2550 |
c ----------------------------------------------- |
| 2551 |
call ctei ( th,q,cloud,cldwat,pl,plk,plke,irun,irise ) |
| 2552 |
|
| 2553 |
return |
| 2554 |
end |
| 2555 |
|
| 2556 |
subroutine ctei ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm ) |
| 2557 |
implicit none |
| 2558 |
integer im,lm |
| 2559 |
real th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm) |
| 2560 |
real plk(im,lm),pl(im,lm),cldfrc(im,lm) |
| 2561 |
integer i,L |
| 2562 |
real getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq |
| 2563 |
real k,krd,kmm,f |
| 2564 |
|
| 2565 |
cp = getcon('CP') |
| 2566 |
alhl = getcon('LATENT HEAT COND') |
| 2567 |
cpoel = cp/alhl |
| 2568 |
elocp = alhl/cp |
| 2569 |
|
| 2570 |
do L=lm,2,-1 |
| 2571 |
do i=1,im |
| 2572 |
dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1) |
| 2573 |
if( dq.eq.0.0 ) dq = 1.0e-20 |
| 2574 |
k = 1.0 + cpoel*plke(i,L)*( th(i,L)-th(i,L-1) ) / dq |
| 2575 |
|
| 2576 |
t = th(i,L)*plk(i,L) |
| 2577 |
p = pl(i,L) |
| 2578 |
call qsat ( t,p,qs,dqsdt,.true. ) |
| 2579 |
|
| 2580 |
krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t ) |
| 2581 |
|
| 2582 |
kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t ) |
| 2583 |
. / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt ) |
| 2584 |
|
| 2585 |
s = ( (k-krd)/(kmm-krd) ) |
| 2586 |
f = 1.0 - min( 1.0, max(0.0,1.0-exp(-s)) ) |
| 2587 |
|
| 2588 |
cldfrc(i,L) = cldfrc(i,L)*f |
| 2589 |
cldwat(i,L) = cldwat(i,L)*f |
| 2590 |
|
| 2591 |
enddo |
| 2592 |
enddo |
| 2593 |
|
| 2594 |
return |
| 2595 |
end |
| 2596 |
|
| 2597 |
subroutine back2grd(gathered,indeces,scattered,irun) |
| 2598 |
implicit none |
| 2599 |
integer i,irun,indeces(irun) |
| 2600 |
real gathered(irun),scattered(irun) |
| 2601 |
real temp(irun) |
| 2602 |
do i = 1,irun |
| 2603 |
temp(indeces(i)) = gathered(i) |
| 2604 |
enddo |
| 2605 |
do i = 1,irun |
| 2606 |
scattered(i) = temp(i) |
| 2607 |
enddo |
| 2608 |
return |
| 2609 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph