File: C:\NOAA\NEMS_11731\src\atmos\phys\rascnvv2.f
1 module module_ras
2 USE MACHINE , ONLY : kind_phys
3 use physcons, grav => con_g, cp => con_cp, alhl => con_hvap &
4 &, alhf => con_hfus, rgas => con_rd, rkap => con_rocp &
5 &, nu => con_FVirt
6 implicit none
7 SAVE
8
9
10
11
12
13
14
15
16
17
18 integer, parameter :: nrcmax=30
19
20 real (kind=kind_phys), parameter :: delt_c=1800.0/3600.0 &
21 &, adjts_d=2.0, adjts_s=0.5
22
23 logical, parameter :: fix_ncld_hr=.true.
24
25 real(kind=kind_phys) ZERO, HALF, ONE, TWO
26 real(kind=kind_phys) FOUR_P2,FOUR
27 real(kind=kind_phys) ONE_M1,ONE_M2,ONE_M5,ONE_M6,ONE_M10
28 PARAMETER (ZERO=0.0, HALF=0.5, ONE=1.0, TWO=2.0)
29 PARAMETER (FOUR_P2=4.E2,FOUR=4.,ONE_M10=1.E-10,ONE_M6=1.E-6 &
30 &, ONE_M5=1.E-5,ONE_M2=1.E-2,ONE_M1=1.E-1)
31
32 real(kind=kind_phys), parameter :: cmb2pa = 100.0
33 real(kind=kind_phys) onebg, gravcon, gravfac, elocp, elfocp, &
34 & rkapi, rkpp1i, zfac, cmpor
35
36 parameter (ONEBG = ONE / GRAV, GRAVCON = cmb2pa * ONEBG &
37 &, GRAVFAC = GRAV / CMB2PA, ELOCP = ALHL / CP &
38 &, ELFOCP = (ALHL+ALHF) / CP &
39 &, RKAPI = ONE / RKAP, RKPP1I = ONE / (ONE+RKAP) &
40 &, CMPOR = CMB2PA / RGAS &
41 &, zfac = 0.28888889E-4 * ONEBG)
42
43
44 logical, parameter :: advcld=.true., advups=.false., advtvd=.true.
45
46 real(kind=kind_phys) RHMAX, qudfac, QUAD_LAM, RHRAM, TESTMB, &
47 & TSTMBI, HCRITD, DD_DP, RKNOB, AFC, EKNOB&
48 &, shalfac,HCRITS, HPERT_FAC
49
50
51 PARAMETER (DD_DP=500.0, RKNOB=1.0, EKNOB=1.0)
52
53
54 PARAMETER (RHMAX=1.0 )
55 PARAMETER (QUAD_LAM=1.0)
56
57 PARAMETER (RHRAM=0.05)
58
59
60 PARAMETER (HCRITD=4000.0)
61 PARAMETER (HCRITS=2000.0)
62
63
64
65
66
67
68
69
70 parameter (qudfac=quad_lam*half, shalfac=3.0)
71
72 parameter (testmb=0.1, tstmbi=one/testmb)
73
74 real(kind=kind_phys) ALMIN1, ALMIN2, ALMAX
75 real(kind=kind_phys) facdt
76
77
78 PARAMETER (ALMIN1=0.00E-6, ALMIN2=0.00E-5, ALMAX=1.0E-2)
79
80
81
82
83 real(kind=kind_phys), parameter :: BLDMAX = 300.0
84
85
86 real(kind=kind_phys) C0, C0I, QI0, QW0, c00, c00i, dlq_fac
87 PARAMETER (QI0=1.0E-5, QW0=1.0E-5)
88
89
90
91 PARAMETER (C00I=1.0E-3)
92
93
94
95 parameter (c00=2.0e-3)
96
97 real(kind=kind_phys) TF, TCR, TCRF, TCL
98
99
100 parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF),TCL=2.0)
101
102
103
104 real(kind=kind_phys) REFP(6), REFR(6), TLAC(8), PLAC(8), TLBPL(7) &
105 &, drdp(5), VTP
106
107 DATA PLAC/100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0/
108 DATA TLAC/ 35.0, 25.0, 20.0, 17.5, 15.0, 12.5, 10.0, 7.5/
109 DATA REFP/500.0, 300.0, 250.0, 200.0, 150.0, 100.0/
110 DATA REFR/ 1.0, 2.0, 3.0, 4.0, 6.0, 8.0/
111
112 real(kind=kind_phys) AC(16), AD(16)
113
114 integer, parameter :: nqrp=500001
115 real(kind=kind_phys) C1XQRP, C2XQRP, TBQRP(NQRP), TBQRA(NQRP) &
116 &, TBQRB(NQRP)
117
118 integer, parameter :: nvtp=10001
119 real(kind=kind_phys) C1XVTP, C2XVTP, TBVTP(NVTP)
120
121 contains
122
123 subroutine set_ras_afc(dt)
124 implicit none
125 real(kind=kind_phys) DT
126
127 = -(1.01097E-4*DT)*(3600./DT)**0.57777778
128 end subroutine set_ras_afc
129
130 subroutine ras_init(levs, me)
131
132 Implicit none
133
134 integer levs
135
136 real(kind=kind_phys) actp, facm, tem, actop, tem1, tem2
137 integer i, l, me
138 PARAMETER (ACTP=1.7, FACM=1.00)
139
140 real(kind=kind_phys) PH(15), A(15)
141
142 DATA PH/150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0 &
143 &, 550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0/
144
145 DATA A/ 1.6851, 1.1686, 0.7663, 0.5255, 0.4100, 0.3677 &
146 &, 0.3151, 0.2216, 0.1521, 0.1082, 0.0750, 0.0664 &
147 &, 0.0553, 0.0445, 0.0633/
148
149 logical first
150 data first/.true./
151
152 if (first) then
153
154 = ACTP*FACM
155 DO L=1,15
156 A(L) = A(L)*FACM
157 ENDDO
158 DO L=2,15
159 TEM = 1.0 / (PH(L) - PH(L-1))
160 AC(L) = (PH(L)*A(L-1) - PH(L-1)*A(L)) * TEM
161 AD(L) = (A(L) - A(L-1)) * TEM
162 ENDDO
163 AC(1) = ACTOP
164 AC(16) = A(15)
165 AD(1) = 0.0
166 AD(16) = 0.0
167
168 CALL SETQRP
169 CALL SETVTP
170
171 do i=1,7
172 tlbpl(i) = (tlac(i)-tlac(i+1)) / (plac(i)-plac(i+1))
173 enddo
174 do i=1,5
175 drdp(i) = (REFR(i+1)-REFR(i)) / (REFP(i+1)-REFP(i))
176 enddo
177
178 = 36.34*SQRT(1.2)* (0.001)**0.1364
179
180 if (me .eq. 0) print *,' NO DOWNDRAFT FOR CLOUD TYPES' &
181 &, ' DETRAINING WITHIN THE BOTTOM ',DD_DP,' hPa LAYERS'
182
183 = .false.
184 endif
185
186 end subroutine ras_init
187 end module module_ras
188
189 module module_rascnv
190
191 USE MACHINE , ONLY : kind_phys
192 implicit none
193 SAVE
194
195 logical REVAP, CUMFRC
196 LOGICAL WRKFUN, CALKBL, CRTFUN, UPDRET, BOTOP, vsmooth
197
198 real(kind=kind_phys), parameter :: frac=0.5, crtmsf=0.0 &
199 &, rhfacs=0.70, rhfacl=0.70 &
200 &, face=5.0, delx=10000.0 &
201 &, ddfac=face*delx*0.001 &
202 &, max_neg_bouy=0.25
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218 PARAMETER ( REVAP=.true., CUMFRC=.true.)
219 PARAMETER (WRKFUN = .FALSE., UPDRET = .FALSE., vsmooth=.false.)
220
221 PARAMETER (CRTFUN = .TRUE., CALKBL = .true., BOTOP=.true.)
222
223
224
225
226
227
228
229
230
231
232
233
234
235 real (kind=kind_phys), parameter :: pgftop=0.0, pgfbot=0.0 &
236
237 , pgfgrad=(pgfbot-pgftop)*0.001
238
239 end module module_rascnv
240
241
242 subroutine rascnv(IM, IX, k, dt, dtf, rannum &
243 &, tin, qin, uin, vin, ccin, trac &
244 &, prsi, prsl, prsik, prslk, phil, phii &
245 &, KPBL, CDRAG, RAINC, kbot, ktop, kcnv &
246 &, DDVEL, FLIPV, facmb, me, garea, lmh, ccwfac&
247 &, nrcm, rhc, ud_mf, dd_mf, det_mf, dlqfac &
248 &, lprnt, ipr, kdt, fscav)
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267 USE MACHINE , ONLY : kind_phys
268 use module_ras, DPD => DD_DP
269 use module_rascnv
270 Implicit none
271
272 LOGICAL FLIPV, lprnt
273
274
275
276 Integer IM, IX, k, ncrnd, me, trac, ipr, nrcm
277 Integer kbot(im), ktop(im), kcnv(im), KPBL(im), lmh(im), kdt
278
279 real(kind=kind_phys) tin(ix,k), qin(ix,k), uin(ix,k) &
280 &, vin(ix,k), prsi(ix,k+1) &
281 &, prsik(ix,k+1), prsl(ix,k), prslk(ix,k+1) &
282 &, phil(ix,k), phii(ix,k+1),ccwfac(im) &
283 &, ccin(ix,k,trac+2) &
284
285 , RAINC(im), CDRAG(im), DDVEL(im) &
286 &, rannum(ix,nrcm),dlqfac &
287 &, ud_mf(im,k), dd_mf(im,k), det_mf(im,k)
288 real(kind=kind_phys) DT, facmb, garea(im), dtf, rhc(im,k) &
289
290
291
292 real(kind=kind_phys), intent(in) :: fscav(trac)
293
294
295
296
297
298
299 real(kind=kind_phys) RAIN, toi(k), qoi(k) &
300 &, TCU(k), QCU(k), PCU(k), clw(k), cli(k) &
301 &, QII(k), QLI(k), PRS(k+1), PSJ(k+1) &
302 &, phi_l(k), phi_h(k+1) &
303 &, wfnc, flx(k+1), FLXD(K+1)
304 real(kind=kind_phys) daylen,pfac,tla,pl,clwmin
305 integer icm,irnd,ib
306
307 PARAMETER (ICM=100, DAYLEN=86400.0, PFAC=1.0/450.0,clwmin=1.0e-10)
308 Integer IC(ICM)
309
310 real(kind=kind_phys), allocatable :: ALFINT(:,:), uvi(:,:)
311 &, trcfac(:,:), rcu(:,:)
312 real(kind=kind_phys) ALFINQ(K), PRSM(K), PSJM(K) &
313 &, alfind(K), rhc_l(k), dtvd(2,4)
314
315 real(kind=kind_phys) CFAC, TEM, dpi, sgc, ccwf, tem1, tem2
316
317 Integer KCR, KFX, NCMX, NC, KTEM, I, L, lm1 &
318 &, ntrc, ia, ll, km1, kp1, ipt, lv, KBL, n &
319 &, lmhij, KRMIN, KRMAX, KFMAX, kblmx
320 real(kind=kind_phys) sgcs(k,im)
321
322 LOGICAL DNDRFT, lprint
323
324
325
326
327 real fscav_(trac+2)
328
329 = 0.0
330 if (trac > 0) then
331 do i=1,trac
332 fscav_(i) = fscav(i)
333 enddo
334 endif
335
336 = k - 1
337 kp1 = k + 1
338
339 = dlqfac
340 tem = 1.0 + dlq_fac
341 c0 = c00 * tem
342 c0i = c00i * tem
343
344 = trac
345 IF (CUMFRC) THEN
346 ntrc = ntrc + 2
347 ENDIF
348 if (ntrc > 0) then
349 if (.not. allocated(trcfac)) allocate (trcfac(k,ntrc))
350 if (.not. allocated(uvi)) allocate (uvi(k,ntrc))
351 if (.not. allocated(rcu)) allocate (rcu(k,ntrc))
352 do n=1, ntrc
353 do l=1,k
354 trcfac(l,n) = 1.0
355 enddo
356 enddo
357 endif
358
359 if (.not. allocated(alfint)) allocate(alfint(k,ntrc+4))
360
361 call set_ras_afc(dt)
362
363 DO IPT=1,IM
364
365 ccwf = 0.5
366 if (ccwfac(ipt) >= 0.0) ccwf = ccwfac(ipt)
367
368
369
370
371
372
373 do l=1,k
374 ud_mf(ipt,l) = 0.0
375 dd_mf(ipt,l) = 0.0
376 det_mf(ipt,l) = 0.0
377 enddo
378
379
380
381
382
383
384 = LMH(ipt)
385 if (flipv) then
386 ll = kp1 - LMHIJ
387 tem = 1.0 / prsi(ipt,ll)
388 else
389 ll = LMHIJ
390 tem = 1.0 / prsi(ipt,ll+1)
391 endif
392 KRMIN = 1
393 KRMAX = km1
394 KFMAX = KRMAX
395 kblmx = 1
396 DO L=1,LMHIJ-1
397 ll = l
398 if (flipv) ll = kp1 -l
399 = prsl(ipt,ll) * tem
400 sgcs(l,ipt) = sgc
401 IF (SGC .LE. 0.050) KRMIN = L
402
403
404 IF (SGC .LE. 0.760) KRMAX = L
405
406 IF (SGC .LE. 0.970) KFMAX = L
407
408 IF (SGC .LE. 0.600) kblmx = L
409
410 ENDDO
411
412
413
414
415 if (fix_ncld_hr) then
416
417 = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/1200) + 0.10001
418
419
420
421
422 = delt_c / dt
423 else
424 NCRND = min(nrcmax, (KRMAX-KRMIN+1))
425 facdt = 1.0 / 3600.0
426 endif
427 NCRND = max(NCRND, 1)
428
429 = MIN(LMHIJ,KRMAX)
430 KTEM = MIN(LMHIJ,KFMAX)
431 KFX = KTEM - KCR
432
433
434
435
436
437
438 IF (KFX .GT. 0) THEN
439 IF (BOTOP) THEN
440 DO NC=1,KFX
441 IC(NC) = KTEM + 1 - NC
442 ENDDO
443 ELSE
444 DO NC=KFX,1,-1
445 IC(NC) = KTEM + 1 - NC
446 ENDDO
447 ENDIF
448 ENDIF
449
450 = KFX + NCRND
451 IF (NCRND .GT. 0) THEN
452 DO I=1,NCRND
453 IRND = (RANNUM(ipt,I)-0.0005)*(KCR-KRMIN+1)
454 IC(KFX+I) = IRND + KRMIN
455 ENDDO
456 ENDIF
457
458
459
460
461
462
463
464
465
466
467
468 = lprnt .and. ipt .eq. ipr
469
470 do l=1,k
471 ll = l
472 if (flipv) ll = kp1 -l
473 (l) = 0.0
474 (l) = 0.0
475
476 (l) = 0.0
477 QLI(l) = 0.0
478
479 (l) = 0.0
480 qcu(l) = 0.0
481 pcu(l) = 0.0
482 flx(l) = 0.0
483 flxd(l) = 0.0
484 rcu(l,1) = 0.0
485 rcu(l,2) = 0.0
486
487 (l) = tin(ipt,ll)
488 qoi(l) = qin(ipt,ll)
489 uvi(l,trac+1) = uin(ipt,ll)
490 uvi(l,trac+2) = vin(ipt,ll)
491
492 if (trac > 0) then
493 do n=1,trac
494 uvi(l,n) = ccin(ipt,ll,n+2)
495 if (abs(uvi(l,n)) < 1.0e-20) uvi(l,n) = 0.0
496 enddo
497 endif
498
499 enddo
500 flx(k+1) = 0.0
501 flxd(k+1) = 0.0
502
503 if (ccin(ipt,1,2) .le. -999.0) then
504 do l=1,k
505 ll = l
506 if (flipv) ll = kp1 -l
507 = ccin(ipt,ll,1) &
508 & * MAX(ZERO, MIN(ONE, (TCR-toi(L))*TCRF))
509 ccin(ipt,ll,2) = ccin(ipt,ll,1) - tem
510 ccin(ipt,ll,1) = tem
511 enddo
512 endif
513 if (advcld) then
514 do l=1,k
515 ll = l
516 if (flipv) ll = kp1 -l
517 (L) = ccin(ipt,ll,1)
518 QLI(L) = ccin(ipt,ll,2)
519 enddo
520 endif
521
522 = KPBL(ipt)
523 if (flipv) KBL = MAX(MIN(k, kp1-KPBL(ipt)), k/2)
524 rain = 0.0
525
526 DO L=1,kp1
527 ll = l
528 if (flipv) ll = kp1 + 1 - l
529 (LL) = prsi(ipt, L) * facmb
530 (LL) = prsik(ipt,L)
531 phi_h(LL) = phii(ipt,L)
532 ENDDO
533
534 DO L=1,k
535 ll = l
536 if (flipv) ll = kp1 - l
537 (LL) = prsl(ipt, L) * facmb
538 (LL) = prslk(ipt,L)
539 phi_l(LL) = phil(ipt,L)
540 rhc_l(LL) = rhc(ipt,L)
541 ENDDO
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560 if (advups) then
561 (:,:) = 1.0
562 elseif (advtvd) then
563 (:,:) = 1.0
564 l = krmin
565 lm1 = l - 1
566 dtvd(1,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1) &
567 & + alhl*(qoi(l)-qoi(lm1))
568 dtvd(1,2) = qoi(l) - qoi(lm1)
569 dtvd(1,3) = qli(l) - qli(lm1)
570 dtvd(1,4) = qii(l) - qii(lm1)
571 do l=krmin+1,k
572 lm1 = l - 1
573
574
575
576
577 (2,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1) &
578 & + alhl*(qoi(l)-qoi(lm1))
579
580
581
582 if (abs(dtvd(2,1)) > 1.0e-10) then
583 tem1 = dtvd(1,1) / dtvd(2,1)
584 tem2 = abs(tem1)
585 alfint(l,1) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)
586 endif
587
588
589
590 (1,1) = dtvd(2,1)
591
592 (2,2) = qoi(l) - qoi(lm1)
593
594
595
596 if (abs(dtvd(2,2)) > 1.0e-10) then
597 tem1 = dtvd(1,2) / dtvd(2,2)
598 tem2 = abs(tem1)
599 alfint(l,2) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)
600 endif
601 dtvd(1,2) = dtvd(2,2)
602
603 (2,3) = qli(l) - qli(lm1)
604
605
606
607 if (abs(dtvd(2,3)) > 1.0e-10) then
608 tem1 = dtvd(1,3) / dtvd(2,3)
609 tem2 = abs(tem1)
610 alfint(l,3) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)
611 endif
612 dtvd(1,3) = dtvd(2,3)
613
614 (2,4) = qii(l) - qii(lm1)
615
616
617
618 if (abs(dtvd(2,4)) > 1.0e-10) then
619 tem1 = dtvd(1,4) / dtvd(2,4)
620 tem2 = abs(tem1)
621 alfint(l,4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)
622 endif
623 dtvd(1,4) = dtvd(2,4)
624 enddo
625
626 if (ntrc > 0) then
627 do n=1,ntrc
628 l = krmin
629 dtvd(1,1) = uvi(l,n) - uvi(l-1,n)
630 do l=krmin+1,k
631 dtvd(2,1) = uvi(l,n) - uvi(l-1,n)
632
633
634
635 if (abs(dtvd(2,1)) > 1.0e-10) then
636 tem1 = dtvd(1,1) / dtvd(2,1)
637 tem2 = abs(tem1)
638 alfint(l,n+4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)
639 endif
640 dtvd(1,1) = dtvd(2,1)
641 enddo
642 enddo
643 endif
644 else
645 alfint(:,:) = 0.5
646 endif
647 alfind(:) = 0.5
648
649
650
651
652
653 if (CUMFRC) then
654 do l=krmin,k
655 tem = 1.0 - max(pgfbot, min(pgftop, pgftop+pgfgrad*prsm(l)))
656 trcfac(l,trac+1) = tem
657 trcfac(l,trac+2) = tem
658 enddo
659 endif
660
661
662
663
664
665
666
667
668
669
670
671 if (calkbl) kbl = k
672 DO NC=1,NCMX
673
674 = IC(NC)
675 if (ib .gt. kbl) cycle
676
677
678
679
680 = DPD > 0.0
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758 = 0.0
759 do L=IB,K+1
760 FLX(L) = 0.0
761 FLXD(L)= 0.0
762 enddo
763
764
765
766
767
768
769
770
771
772
773 = -10.0
774
775 CALL CLOUD(lmhij, IB, ntrc, kblmx &
776 &, FRAC, MAX_NEG_BOUY, vsmooth &
777 &, REVAP, WRKFUN, CALKBL, CRTFUN, DNDRFT, lprint &
778 &, DT, KDT, TLA, DPD &
779 &, ALFINT, rhfacl, rhfacs, garea(ipt) &
780 &, ccwf, CDRAG(ipt), trcfac &
781 &, alfind, rhc_l, phi_l, phi_h, PRS, PRSM,sgcs(1,ipt) &
782 &, TOI, QOI, UVI, QLI, QII, KBL, DDVEL(ipt) &
783 &, TCU, QCU, RCU, PCU, FLX, FLXD, RAIN, WFNC, fscav_ &
784 & )
785
786
787
788
789
790
791
792
793
794
795
796
797
798 do L=IB,K
799 ll = l
800 if (flipv) ll = kp1 -l
801 (ipt,ll) = ud_mf(ipt,ll) + flx(l+1)
802 dd_mf(ipt,ll) = dd_mf(ipt,ll) + flxd(l+1)
803 enddo
804 ll = ib
805 if (flipv) ll = kp1 - ib
806 det_mf(ipt,ll) = det_mf(ipt,ll) + flx(ib)
807
808
809
810
811
812
813
814 if (.not. advcld) then
815 do l=1,K
816 clw(l ) = clw(l) + QLI(L)
817 cli(l ) = cli(l) + QII(L)
818 QLI(L) = 0.0
819 QII(L) = 0.0
820 enddo
821 endif
822
823 ENDDO
824
825 (ipt) = rain * 0.001
826
827
828
829
830
831
832
833
834 do l=1,k
835 ll = l
836 if (flipv) ll = kp1 - l
837 tin(ipt,ll) = toi(l)
838 (ipt,ll) = qoi(l)
839 (ipt,ll) = uvi(l,trac+1)
840 (ipt,ll) = uvi(l,trac+2)
841 if (trac > 0) then
842 do n=1,trac
843 ccin(ipt,ll,n+2) = uvi(l,n)
844 enddo
845 endif
846 enddo
847 if (advcld) then
848 do l=1,k
849 ll = l
850 if (flipv) ll = kp1 - l
851 ccin(ipt,ll,1) = qii(l)
852 (ipt,ll,2) = qli(l)
853 enddo
854 else
855 do l=1,k
856 ll = l
857 if (flipv) ll = kp1 - l
858 ccin(ipt,ll,1) = ccin(ipt,ll,1) + cli(l)
859 ccin(ipt,ll,2) = ccin(ipt,ll,2) + clw(l)
860 enddo
861 endif
862
863 (ipt) = kp1
864 kbot(ipt) = 0
865
866 kcnv(ipt) = 0
867
868 do l=lmhij-1,1,-1
869 if (sgcs(l,ipt) < 0.85 .and. tcu(l) .ne. 0.0) then
870
871 (ipt) = 1
872 endif
873
874 if (clw(l)+cli(l) .gt. 0.0 .OR. &
875 & qli(l)+qii(l) .gt. clwmin) ktop(ipt) = l
876 enddo
877 do l=1,km1
878 if (clw(l)+cli(l) .gt. 0.0 .OR. &
879 & qli(l)+qii(l) .gt. clwmin) kbot(ipt) = l
880 enddo
881 if (flipv) then
882 ktop(ipt) = kp1 - ktop(ipt)
883 kbot(ipt) = kp1 - kbot(ipt)
884 endif
885
886
887
888
889
890
891
892
893 (ipt) = DDVEL(ipt) * DDFAC * GRAV / (prs(K+1)-prs(k))
894
895 ENDDO
896
897 deallocate (alfint,uvi,trcfac,rcu)
898
899 RETURN
900 END
901 SUBROUTINE CRTWRK(PL, CCWF, ACR)
902 USE MACHINE , ONLY : kind_phys
903 use module_ras , only : ac, ad
904 Implicit none
905
906 real(kind=kind_phys) PL, CCWF, ACR
907 INTEGER IWK
908
909 = PL * 0.02 - 0.999999999
910 IWK = MAX(1, MIN(IWK,16))
911 ACR = (AC(IWK) + PL * AD(IWK)) * CCWF
912
913 RETURN
914 END
915 SUBROUTINE CLOUD( &
916 & K, KD, NTRC, KBLMX &
917 &, FRACBL, MAX_NEG_BOUY, vsmooth &
918 &, REVAP, WRKFUN, CALKBL, CRTFUN, DNDRFT, lprnt &
919 &, DT, KDT, TLA, DPD &
920 &, ALFINT, RHFACL, RHFACS, garea, ccwf, cd, trcfac &
921 &, alfind, rhc_ls, phil, phih, prs, prsm, sgcs &
922 &, TOI, QOI, ROI, QLI, QII, KPBL, DSFC &
923 &, TCU, QCU, RCU, PCU, FLX, FLXD, CUP, WFNC,fscav_ &
924 & )
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978 USE MACHINE , ONLY : kind_phys
979 use module_ras
980 IMPLICIT NONE
981
982
983
984
985 LOGICAL REVAP, DNDRFT, WRKFUN, CALKBL, CRTFUN, CALCUP
986 logical vsmooth, lprnt
987 INTEGER K, KD, NTRC, kblmx
988
989
990 real(kind=kind_phys) TOI(K), QOI(K ), PRS(K+1), PRSM(K) &
991 &, QLI(K), QII(K), PHIH(K+1), PHIL(K) &
992 &, ROI(K,NTRC), SGCS(K)
993 real(kind=kind_phys) CD, UFN, DSFC
994 INTEGER KPBL, KBL, KB1, kdt
995
996 real(kind=kind_phys) FRACBL, MAX_NEG_BOUY, &
997 & ALFINT(K,NTRC+4), RHFACL, RHFACS, garea, ccwf
998 real(kind=kind_phys) DPD, alfind(k), rhc_ls(k)
999 real(kind=kind_phys) trcfac(k,NTRC)
1000
1001
1002
1003 real(kind=kind_phys) TCU(K), QCU(K), RCU(K,NTRC) &
1004 &, TCD(K), QCD(K), PCU(K) &
1005 &, FLX(K+1), FLXD(K+1), CUP
1006
1007
1008
1009 real(kind=kind_phys) HOL(KD:K), QOL(KD:K), GAF(KD:K+1) &
1010 &, HST(KD:K), QST(KD:K), TOL(KD:K) &
1011 &, GMH(KD:K), GMS(KD:K+1), GAM(KD:K+1) &
1012 &, AKT(KD:K), AKC(KD:K), BKC(KD:K) &
1013 &, LTL(KD:K), RNN(KD:K), FCO(KD:K) &
1014 &, PRI(KD:K) &
1015 &, QIL(KD:K), QLL(KD:K) &
1016 &, ZET(KD:K), XI(KD:K), RNS(KD:K) &
1017 &, Q0U(KD:K), Q0D(KD:K), vtf(KD:K) &
1018 &, DLB(KD:K+1),DLT(KD:K+1), ETA(KD:K+1) &
1019 &, PRL(KD:K+1) &
1020 &, CIL(KD:K), CLL(KD:K), ETAI(KD:K) &
1021 &, dlq(kd:k)
1022
1023 real(kind=kind_phys) ALM, DET, HCC, CLP &
1024 &, HSU, HSD, QTL, QTV &
1025 &, AKM, WFN, HOS, QOS &
1026 &, AMB, TX1, TX2, TX3 &
1027 &, TX4, TX5, QIS, QLS &
1028 &, HBL, QBL, RBL(NTRC) &
1029 &, QLB, QIB, PRIS &
1030 &, WFNC, TX6, ACR &
1031 &, TX7, TX8, TX9, RHC &
1032 &, hstkd, qstkd, ltlkd, q0ukd, q0dkd, dlbkd &
1033 &, qtp, qw00, qi00, qrbkd &
1034 &, hstold, rel_fac, prism
1035
1036 real(kind=kind_phys) fscav_(ntrc)
1037
1038 real(kind=kind_phys) wrk1(kd:k), wrk2(kd:k)
1039
1040
1041 LOGICAL ep_wfn, cnvflg
1042
1043 LOGICAL LOWEST, SKPDD
1044
1045 real(kind=kind_phys) TL, PL, QL, QS, DQS, ST1, SGN, TAU, &
1046 & QTVP, HB, QB, TB, QQQ, &
1047 & HCCP, DS, DH, AMBMAX, X00, EPP, QTLP, &
1048 & DPI, DPHIB, DPHIT, DEL_ETA, DETP, &
1049 & TEM, TEM1, TEM2, TEM3, TEM4, &
1050 & ST2, ST3, ST4, ST5, &
1051 & ERRMIN, ERRMI2, ERRH, ERRW, ERRE, TEM5, &
1052 & TEM6, HBD, QBD, st1s, shal_fac, hmax, hmin, &
1053 & dhdpmn, dhdp(kd:k)
1054 parameter (ERRMIN=0.0001, ERRMI2=0.1*ERRMIN)
1055 INTEGER I, L, N, KD1, II &
1056 &, KP1, IT, KM1, KTEM, KK, KK1, LM1, LL, LP1, kbls, kmxh
1057 &, kblh, kblm, kblpmn, kmax, kmaxm1, kmaxp1, klcl, kmin, kmxb
1058
1059 real avt, avq, avr, avh
1060
1061
1062
1063 real(kind=kind_phys), parameter :: rainmin=1.0e-8
1064 real(kind=kind_phys), parameter :: oneopt9=1.0/0.09
1065 real(kind=kind_phys), parameter :: oneopt4=1.0/0.04
1066
1067 real(kind=kind_phys) CLFRAC, DT, clf, clvfr
1068
1069 real(kind=kind_phys) ACTEVAP,AREARAT,DELTAQ,MASS,MASSINV,POTEVAP &
1070 &, TEQ,QSTEQ,DQDT,QEQ
1071
1072
1073
1074 real(kind=kind_phys) TLA, GMF
1075
1076 real(kind=kind_phys) BUY(KD:K+1), QRB(KD:K), QRT(KD:K) &
1077 &, ETD(KD:K+1), HOD(KD:K+1), QOD(KD:K+1) &
1078 &, GHD(KD:K), GSD(KD:K), EVP(KD:K) &
1079 &, ETZ(KD:K), CLDFR(KD:K), ETZI(KD:K-1) &
1080 &, TRAIN, DOF, CLDFRD &
1081 &, FAC, RSUM1, RSUM2, RSUM3, dpneg, hcrit
1082 INTEGER IDH, lcon
1083 LOGICAL DDFT, UPDRET
1084
1085
1086
1087 real delzkm
1088 real fnoscav
1089
1090
1091
1092 do l=1,K
1093 tcd(L) = 0.0
1094 qcd(L) = 0.0
1095 enddo
1096
1097 = K + 1
1098 KM1 = K - 1
1099 KD1 = KD + 1
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111 = 0.0
1112 DOF = 0.0
1113 PRL(KP1) = PRS(KP1)
1114
1115 DO L=KD,K
1116 RNN(L) = 0.0
1117 ZET(L) = 0.0
1118 XI(L) = 0.0
1119
1120 (L) = TOI(L)
1121 QOL(L) = QOI(L)
1122 PRL(L) = PRS(L)
1123 BUY(L) = 0.0
1124 CLL(L) = QLI(L)
1125 CIL(L) = QII(L)
1126 ENDDO
1127
1128 if (vsmooth) then
1129 do l=kd,k
1130 wrk1(l) = tol(l)
1131 wrk2(l) = qol(l)
1132 enddo
1133 do l=kd1,km1
1134 tol(l) = 0.25*wrk1(l-1) + 0.5*wrk1(l) + 0.25*wrk1(l+1)
1135 qol(l) = 0.25*wrk2(l-1) + 0.5*wrk2(l) + 0.25*wrk2(l+1)
1136 enddo
1137 endif
1138
1139 DO L=KD, K
1140 DPI = ONE / (PRL(L+1) - PRL(L))
1141 PRI(L) = GRAVFAC * DPI
1142
1143 = PRSM(L)
1144 TL = TOL(L)
1145
1146 AKT(L) = (PRL(L+1) - PL) * DPI
1147
1148 CALL QSATCN(TL, PL, QS, DQS)
1149
1150
1151 (L) = QS
1152 GAM(L) = DQS * ELOCP
1153 ST1 = ONE + GAM(L)
1154 GAF(L) = (ONE/ALHL) * GAM(L)/ST1
1155
1156 QL = MAX(MIN(QS*RHMAX,QOL(L)), ONE_M10)
1157 QOL(L) = QL
1158
1159 TEM = CP * TL
1160 LTL(L) = TEM * ST1 / (ONE+NU*(QST(L)+TL*DQS))
1161 vtf(L) = 1.0 + NU * QL
1162 ETA(L) = ONE / (LTL(L) * VTF(L))
1163
1164 HOL(L) = TEM + QL * ALHL
1165 HST(L) = TEM + QS * ALHL
1166
1167 ENDDO
1168
1169 (K+1) = ZERO
1170 GMS(K) = ZERO
1171
1172 (KD) = HALF
1173 GMS(KD) = ZERO
1174
1175 = ZERO
1176
1177 (K+1) = GAM(K)
1178 GAF(K+1) = GAF(K)
1179
1180 DO L=K,KD1,-1
1181 DPHIB = PHIL(L) - PHIH(L+1)
1182 DPHIT = PHIH(L) - PHIL(L)
1183
1184 (L) = DPHIB * ETA(L)
1185 DLT(L) = DPHIT * ETA(L)
1186
1187 (L) = DPHIB
1188 QRT(L) = DPHIT
1189
1190 (L) = ETA(L+1) + DPHIB
1191
1192 HOL(L) = HOL(L) + ETA(L)
1193 hstold = hst(l)
1194 HST(L) = HST(L) + ETA(L)
1195
1196 (L) = ETA(L) + DPHIT
1197 ENDDO
1198
1199
1200
1201 = KD
1202
1203 DPHIB = PHIL(L) - PHIH(L+1)
1204
1205 (L) = DPHIB * ETA(L)
1206
1207 (L) = DPHIB
1208 QRT(L) = DPHIB
1209
1210 (L) = ETA(L+1) + DPHIB
1211
1212 HOL(L) = HOL(L) + ETA(L)
1213 HST(L) = HST(L) + ETA(L)
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231 = hcritd
1232 if (sgcs(kd) > 0.65) hcrit = hcrits
1233 IF (CALKBL) THEN
1234 KTEM = MAX(KD+1, KBLMX)
1235 hmin = hol(k)
1236 kmin = k
1237 do l=km1,kd,-1
1238 if (hmin > hol(l)) then
1239 hmin = hol(l)
1240 kmin = l
1241 endif
1242 enddo
1243 if (kmin == k) return
1244 hmax = hol(k)
1245 kmax = k
1246 do l=km1,ktem,-1
1247 if (hmax < hol(l)) then
1248 hmax = hol(l)
1249 kmax = l
1250 endif
1251 enddo
1252 kmxb = kmax
1253 if (kmax < kmin) then
1254 kmax = k
1255 kmxb = k
1256 hmax = hol(kmax)
1257 elseif (kmax < k) then
1258 do l=kmax+1,k
1259 if (abs(hol(kmax)-hol(l)) > 0.5 * hcrit) then
1260 kmxb = l - 1
1261 exit
1262 endif
1263 enddo
1264 endif
1265 kmaxm1 = kmax - 1
1266 kmaxp1 = kmax + 1
1267 kblpmn = kmax
1268
1269 (kmax:k) = 0.0
1270 dhdpmn = dhdp(kmax)
1271 do l=kmaxm1,ktem,-1
1272 dhdp(l) = (HOL(L)-HOL(L+1)) / (PRL(L+2)-PRL(L))
1273 if (dhdp(l) < dhdpmn) then
1274 dhdpmn = dhdp(l)
1275 kblpmn = l + 1
1276 elseif (dhdp(l) > 0.0 .and. l <= kmin) then
1277 exit
1278 endif
1279 enddo
1280 kbl = kmax
1281 if (kblpmn < kmax) then
1282 do l=kblpmn,kmaxm1
1283 if (hmax-hol(l) < 0.5*hcrit) then
1284 kbl = l
1285 exit
1286 endif
1287 enddo
1288 endif
1289
1290
1291
1292 = kd1
1293 if (kmax > kd1) then
1294 do l=kmaxm1,kd1,-1
1295 if (hmax > hst(l)) then
1296 klcl = l+1
1297 exit
1298 endif
1299 enddo
1300 endif
1301
1302
1303
1304
1305
1306 if (kmax < kmxb) then
1307 kmax = max(kd1, min(kmxb,k))
1308 kmaxm1 = kmax - 1
1309 kmaxp1 = kmax + 1
1310 endif
1311
1312
1313
1314
1315 = max(kbl,kd1)
1316 kbl = max(klcl,kd1)
1317 if (prl(kmaxp1) - prl(ii) > 50.0 .and. ii > kbl) kbl = ii
1318
1319
1320
1321 if (kbl .ne. ii) then
1322 if (PRL(kmaxp1)-PRL(KBL) > bldmax) kbl = max(kbl,ii)
1323 endif
1324 if (kbl < ii) then
1325 if (hol(ii)-hol(ii-1) > 0.5*hcrit) kbl = ii
1326 endif
1327
1328
1329 if (prl(kbl) - prl(klcl) > 250.0 ) return
1330
1331 = min(kmax, MAX(KBL,KBLMX))
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351 = KBL
1352
1353
1354
1355
1356
1357 ELSE
1358 KBL = KPBL
1359
1360 ENDIF
1361
1362
1363
1364
1365 = min(kmax,MAX(KBL,KD+2))
1366 KB1 = KBL - 1
1367
1368
1369
1370 if(PRL(Kmaxp1)-PRL(KBL) .gt. bldmax .or. kb1 .le. kd) then
1371 return
1372 endif
1373
1374
1375
1376
1377 = ONE / (PRL(K+1)-PRL(KBL))
1378 PRISM = ONE / (PRL(Kmaxp1)-PRL(KBL))
1379 TX1 = ETA(KBL)
1380
1381 (KBL) = 0.0
1382 XI(KBL) = 0.0
1383 ZET(KBL) = 0.0
1384
1385 = 1.0
1386
1387 if (prl(kbl)-prl(kd) < 350.0 .and. kmax == k) shal_fac = shalfac
1388 DO L=Kmax,KD,-1
1389 IF (L >= KBL) THEN
1390 ETA(L) = (PRL(Kmaxp1)-PRL(L)) * PRISM
1391 ELSE
1392 ZET(L) = (ETA(L) - TX1) * ONEBG
1393 XI(L) = ZET(L) * ZET(L) * (QUDFAC*shal_fac)
1394 ETA(L) = ZET(L) - ZET(L+1)
1395 GMS(L) = XI(L) - XI(L+1)
1396 ENDIF
1397
1398 ENDDO
1399 if (kmax < k) then
1400 do l=kmaxp1,kp1
1401 eta(l) = 0.0
1402 enddo
1403 endif
1404
1405 = HOL(Kmax) * ETA(Kmax)
1406 QBL = QOL(Kmax) * ETA(Kmax)
1407 QLB = CLL(Kmax) * ETA(Kmax)
1408 QIB = CIL(Kmax) * ETA(Kmax)
1409 TX1 = QST(Kmax) * ETA(Kmax)
1410
1411 DO L=Kmaxm1,KBL,-1
1412 TEM = ETA(L) - ETA(L+1)
1413 HBL = HBL + HOL(L) * TEM
1414 QBL = QBL + QOL(L) * TEM
1415 QLB = QLB + CLL(L) * TEM
1416 QIB = QIB + CIL(L) * TEM
1417 TX1 = TX1 + QST(L) * TEM
1418 ENDDO
1419
1420
1421
1422
1423
1424
1425
1426
1427 = HOL(KD)
1428 IDH = KD1
1429 DO L=KD1,KB1
1430 IF (HOL(L) .LT. TX2) THEN
1431 TX2 = HOL(L)
1432 IDH = L
1433 ENDIF
1434 ENDDO
1435 IDH = 1
1436 IDH = MAX(KD1, IDH)
1437
1438 = HBL - HOL(KD)
1439 TEM = HBL - HST(KD1) - LTL(KD1) * NU *(QOL(KD1)-QST(KD1))
1440 LOWEST = KD .EQ. KB1
1441
1442 lcon = kd
1443 do l=kb1,kd1,-1
1444 if (hbl >= hst(l)) then
1445 lcon = l
1446 exit
1447 endif
1448 enddo
1449
1450 if (lcon == kd .or. kbl <= kd .or. prl(kbl)-prsm(lcon) > 150.0) &
1451 & return
1452
1453 = RHFACS - QBL / TX1
1454
1455 = (TEM .GT. ZERO .OR. (LOWEST .AND. TEM1 .GE. ZERO)) &
1456 & .AND. (TX1 .LT. RHRAM)
1457
1458
1459
1460
1461
1462
1463
1464
1465 IF (.NOT. cnvflg) RETURN
1466
1467 = MAX(ZERO, MIN(ONE, EXP(-20.0*TX1) ))
1468
1469 if (ntrc > 0) then
1470 DO N=1,NTRC
1471 RBL(N) = ROI(Kmax,N) * ETA(Kmax)
1472 ENDDO
1473 DO N=1,NTRC
1474 DO L=KmaxM1,KBL,-1
1475 RBL(N) = RBL(N) + ROI(L,N)*(ETA(L)-ETA(L+1))
1476 ENDDO
1477 ENDDO
1478 endif
1479
1480 = 0.0
1481 TX5 = 0.0
1482
1483 = QST(KBL) - GAF(KBL) * HST(KBL)
1484 DO L=KBL,K
1485 QIL(L) = MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(L))*TCRF))
1486 ENDDO
1487
1488 DO L=KB1,KD1,-1
1489 TEM = QST(L) - GAF(L) * HST(L)
1490 TEM1 = (TX3 + TEM) * 0.5
1491 ST2 = (GAF(L)+GAF(L+1)) * 0.5
1492
1493 (L+1) = TEM1 + ST2 * HBL
1494
1495
1496
1497
1498 (L+1) = ZET(L+1) * TEM1 + ST2 * TX4
1499 GMH(L+1) = XI(L+1) * TEM1 + ST2 * TX5
1500
1501 = TEM
1502 TX4 = TX4 + ETA(L) * HOL(L)
1503 TX5 = TX5 + GMS(L) * HOL(L)
1504
1505 (L) = MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(L))*TCRF))
1506 QLL(L+1) = (0.5*ALHF) * ST2 * (QIL(L)+QIL(L+1)) + ONE
1507 ENDDO
1508
1509
1510
1511 = KD
1512
1513 = QST(L) - GAF(L) * HST(L)
1514 TEM1 = (TX3 + TEM) * 0.5
1515 ST2 = (GAF(L)+GAF(L+1)) * 0.5
1516
1517 (L+1) = TEM1 + ST2 * HBL
1518 RNN(L+1) = ZET(L+1) * TEM1 + ST2 * TX4
1519 GMH(L+1) = XI(L+1) * TEM1 + ST2 * TX5
1520
1521 (L) = TEM + GAF(L) * HBL
1522 RNN(L) = TEM * ZET(L) + (TX4 + ETA(L)*HOL(L)) * GAF(L)
1523 GMH(L) = TEM * XI(L) + (TX5 + GMS(L)*HOL(L)) * GAF(L)
1524
1525
1526
1527 (KBL) = QBL
1528 RNN(KBL) = 0.0
1529 GMH(KBL) = 0.0
1530
1531 (KD) = MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(KD))*TCRF))
1532 QLL(KD1) = (0.5*ALHF) * ST2 * (QIL(KD) + QIL(KD1)) + ONE
1533 QLL(KD ) = ALHF * GAF(KD) * QIL(KD) + ONE
1534
1535
1536
1537
1538
1539
1540
1541 = qil(kd)
1542 st2 = c0i * st1
1543 tem = c0 * (1.0-st1)
1544 tem2 = st2*qi0 + tem*qw0
1545
1546 DO L=KD,KB1
1547 tx2 = akt(l) * eta(l)
1548 tx1 = tx2 * tem2
1549 q0u(l) = tx1
1550 FCO(L) = FCO(L+1) - FCO(L) + tx1
1551 RNN(L) = RNN(L+1) - RNN(L) &
1552 & + ETA(L)*(QOL(L)+CLL(L)+CIL(L)) + tx1*zet(l)
1553 GMH(L) = GMH(L+1) - GMH(L) &
1554 & + GMS(L)*(QOL(L)+CLL(L)+CIL(L)) + tx1*xi(l)
1555
1556 = (1.0-akt(l)) * eta(l)
1557
1558
1559
1560
1561 (L) = QLL(L) + (st2 + tem) * tx2
1562
1563
1564
1565 (L) = 1.0 / AKT(L)
1566
1567 = 0.5 * (qil(l)+qil(l+1))
1568 st2 = c0i * st1
1569 tem = c0 * (1.0-st1)
1570 tem2 = st2*qi0 + tem*qw0
1571
1572 (L) = QLL(L+1) - (st2 + tem) * tem1
1573
1574 = tem1*tem2
1575 q0d(l) = tx1
1576 FCO(L) = FCO(L) + tx1
1577 RNN(L) = RNN(L) + tx1*zet(l+1)
1578 GMH(L) = GMH(L) + tx1*xi(l+1)
1579 ENDDO
1580
1581
1582
1583
1584 = qw0
1585 qi00 = qi0
1586 ii = 0
1587 777 continue
1588
1589
1590
1591 = .false.
1592 RNN(KBL) = 0.0
1593 TX3 = bkc(kb1) * (QIB + QLB)
1594 TX4 = 0.0
1595 TX5 = 0.0
1596 DO L=KB1,KD1,-1
1597 TEM = BKC(L-1) * AKC(L)
1598
1599
1600 = (TX3 + FCO(L)) * TEM
1601 TX4 = (TX4 + RNN(L)) * TEM
1602 TX5 = (TX5 + GMH(L)) * TEM
1603 ENDDO
1604 IF (KD .LT. KB1) THEN
1605 HSD = HST(KD1) + LTL(KD1) * NU *(QOL(KD1)-QST(KD1))
1606 ELSE
1607 HSD = HBL
1608 ENDIF
1609
1610
1611
1612 = (TX3 + FCO(KD)) * AKC(KD)
1613 TX4 = (TX4 + RNN(KD)) * AKC(KD)
1614 TX5 = (TX5 + GMH(KD)) * AKC(KD)
1615 ALM = ALHF*QIL(KD) - LTL(KD) * VTF(KD)
1616
1617 = HST(KD) + LTL(KD) * NU * (QOL(KD)-QST(KD))
1618
1619
1620
1621
1622
1623
1624 = ALM * TX4
1625 TX2 = ALM * TX5
1626
1627 DO L=KD,KB1
1628 TAU = HOL(L) - HSU
1629 TX1 = TX1 + TAU * ETA(L)
1630 TX2 = TX2 + TAU * GMS(L)
1631 ENDDO
1632
1633
1634
1635
1636
1637 = HSU - ALM * TX3
1638
1639 = ZERO
1640 ALM = -100.0
1641 HOS = HOL(KD)
1642 QOS = QOL(KD)
1643 QIS = CIL(KD)
1644 QLS = CLL(KD)
1645 cnvflg = HBL .GT. HSU .and. abs(tx1) .gt. 1.0e-4
1646
1647
1648
1649
1650
1651
1652
1653
1654 = HALF*(HSU + HSD)
1655 IF (cnvflg) THEN
1656
1657
1658
1659
1660
1661 = 1.0
1662 st2 = hbl - hsu
1663
1664
1665
1666 if (tx2 .eq. 0.0) then
1667 alm = - st2 / tx1
1668 if (alm .gt. almax) alm = -100.0
1669 else
1670 x00 = tx2 + tx2
1671 epp = tx1 * tx1 - (x00+x00)*st2
1672 if (epp .gt. 0.0) then
1673 x00 = 1.0 / x00
1674 tem = sqrt(epp)
1675 tem1 = (-tx1-tem)*x00
1676 tem2 = (-tx1+tem)*x00
1677 if (tem1 .gt. almax) tem1 = -100.0
1678 if (tem2 .gt. almax) tem2 = -100.0
1679 alm = max(tem1,tem2)
1680
1681
1682
1683
1684 endif
1685 endif
1686
1687
1688
1689
1690
1691
1692
1693
1694 ELSEIF ( (HBL .LE. HSU) .AND. &
1695 & (HBL .GT. ST1 ) ) THEN
1696 ALM = ZERO
1697
1698 ENDIF
1699
1700 = .TRUE.
1701 IF (ALMIN1 .GT. 0.0) THEN
1702 IF (ALM .GE. ALMIN1) cnvflg = .FALSE.
1703 ELSE
1704 LOWEST = KD .EQ. KB1
1705 IF ( (ALM .GT. ZERO) .OR. &
1706 & (.NOT. LOWEST .AND. ALM .EQ. ZERO) ) cnvflg = .FALSE.
1707 ENDIF
1708
1709
1710
1711 IF (cnvflg) THEN
1712 IF (ii .gt. 0 .or. (qw00 .eq. 0.0 .and. qi00 .eq. 0.0)) RETURN
1713 CLP = 1.0
1714 ep_wfn = .true.
1715 GO TO 888
1716 ENDIF
1717
1718
1719
1720
1721 = ONE
1722 IF(CLP.GT.ZERO .AND. CLP.LT.ONE) THEN
1723 ST1 = HALF*(ONE+CLP)
1724 ST2 = ONE - ST1
1725 st1s = st1
1726 hstkd = hst(kd)
1727 qstkd = qst(kd)
1728 ltlkd = ltl(kd)
1729 q0ukd = q0u(kd)
1730 q0dkd = q0d(kd)
1731 dlbkd = dlb(kd)
1732 qrbkd = qrb(kd)
1733
1734 (KD) = HST(KD)*ST1 + HST(KD1)*ST2
1735 HOS = HOL(KD)*ST1 + HOL(KD1)*ST2
1736 QST(KD) = QST(KD)*ST1 + QST(KD1)*ST2
1737 QOS = QOL(KD)*ST1 + QOL(KD1)*ST2
1738 QLS = CLL(KD)*ST1 + CLL(KD1)*ST2
1739 QIS = CIL(KD)*ST1 + CIL(KD1)*ST2
1740 LTL(KD) = LTL(KD)*ST1 + LTL(KD1)*ST2
1741
1742 (KD) = DLB(KD)*CLP
1743 qrb(KD) = qrb(KD)*CLP
1744 ETA(KD) = ETA(KD)*CLP
1745 GMS(KD) = GMS(KD)*CLP
1746 Q0U(KD) = Q0U(KD)*CLP
1747 Q0D(KD) = Q0D(KD)*CLP
1748 ENDIF
1749
1750
1751
1752
1753
1754
1755 = 0.0
1756 TEM = PRL(KD1) - (PRL(KD1)-PRL(KD)) * CLP * HALF
1757 tx1 = PRL(KBL) - TEM
1758 tx2 = min(900.0,max(tx1,100.0))
1759 tem1 = log(tx2*0.01) / log(10.0)
1760 if ( kdt == 1 ) then
1761 rel_fac = (dt * facdt) / (tem1*12.0 + (1-tem1)*3.0)
1762 else
1763 rel_fac = (dt * facdt) / (tem1*adjts_d + (1-tem1)*adjts_s)
1764 endif
1765
1766 = max(zero, min(one,rel_fac))
1767
1768 IF (CRTFUN) THEN
1769 CALL CRTWRK(TEM, CCWF, ST1)
1770 ACR = TX1 * ST1
1771 ENDIF
1772
1773
1774
1775
1776
1777
1778
1779
1780 DO L=KB1,KD,-1
1781 ETA(L) = ETA(L+1) + ALM * (ETA(L) + ALM * GMS(L))
1782 ENDDO
1783 DO L=KD,KBL
1784 ETAI(L) = 1.0 / ETA(L)
1785 ENDDO
1786
1787
1788
1789
1790
1791 = ZERO
1792 AKM = ZERO
1793 DET = ZERO
1794 HCC = HBL
1795 cnvflg = .FALSE.
1796 QTL = QST(KB1) - GAF(KB1)*HST(KB1)
1797 TX1 = HBL
1798
1799 = qbl
1800 det = qlb + qib
1801
1802 = 0.0
1803 dpneg = 0.0
1804
1805 DO L=KB1,KD1,-1
1806 DEL_ETA = ETA(L) - ETA(L+1)
1807 HCCP = HCC + DEL_ETA*HOL(L)
1808
1809 = QST(L-1) - GAF(L-1)*HST(L-1)
1810 QTVP = 0.5 * ((QTLP+QTL)*ETA(L) &
1811 & + (GAF(L)+GAF(L-1))*HCCP)
1812 ST1 = ETA(L)*Q0U(L) + ETA(L+1)*Q0D(L)
1813 DETP = (BKC(L)*DET - (QTVP-QTV) &
1814 & + DEL_ETA*(QOL(L)+CLL(L)+CIL(L)) + ST1) * AKC(L)
1815
1816
1817
1818
1819
1820
1821
1822 = AKT(L) - QLL(L)
1823 TEM2 = QLL(L+1) - BKC(L)
1824 RNS(L) = TEM1*DETP + TEM2*DET - ST1
1825
1826 qtp = 0.5 * (qil(L)+qil(L-1))
1827 tem2 = min(qtp*(detp-eta(l)*qw00), &
1828 & (1.0-qtp)*(detp-eta(l)*qi00))
1829 st1 = min(tx2,tem2)
1830 tx2 = tem2
1831
1832 IF (rns(l) .lt. zero .or. st1 .lt. zero) ep_wfn = .TRUE.
1833 IF (DETP <= ZERO) cnvflg = .TRUE.
1834
1835 ST1 = HST(L) - LTL(L)*NU*(QST(L)-QOL(L))
1836
1837
1838 TEM2 = HCCP + DETP * QTP * ALHF
1839
1840
1841
1842
1843
1844
1845
1846 = LTL(L) * VTF(L)
1847 TEM5 = CLL(L) + CIL(L)
1848 TEM3 = (TX1 - ETA(L+1)*ST1 - ST2*(DET-TEM5*eta(l+1))) * DLB(L)
1849 TEM4 = (TEM2 - ETA(L )*ST1 - ST2*(DETP-TEM5*eta(l))) * DLT(L)
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861 = TEM3 + TEM4
1862
1863
1864
1865
1866 = WFN + ST1
1867 AKM = AKM - min(ST1,ZERO)
1868
1869
1870
1871 if (st1 .lt. zero .and. wfn .lt. zero) then
1872 dpneg = dpneg + prl(l+1) - prl(l)
1873 endif
1874
1875 BUY(L) = 0.5 * (tem3/(eta(l+1)*qrb(l)) + tem4/(eta(l)*qrt(l)))
1876
1877 = HCCP
1878 DET = DETP
1879 QTL = QTLP
1880 QTV = QTVP
1881 TX1 = TEM2
1882
1883 ENDDO
1884
1885 DEL_ETA = ETA(KD) - ETA(KD1)
1886 HCCP = HCC + DEL_ETA*HOS
1887
1888 = QST(KD) - GAF(KD)*HST(KD)
1889 QTVP = QTLP*ETA(KD) + GAF(KD)*HCCP
1890 ST1 = ETA(KD)*Q0U(KD) + ETA(KD1)*Q0D(KD)
1891 DETP = (BKC(KD)*DET - (QTVP-QTV) &
1892 & + DEL_ETA*(QOS+QLS+QIS) + ST1) * AKC(KD)
1893
1894 = AKT(KD) - QLL(KD)
1895 TEM2 = QLL(KD1) - BKC(KD)
1896 RNS(KD) = TEM1*DETP + TEM2*DET - ST1
1897
1898 IF (rns(kd) .lt. zero) ep_wfn = .TRUE.
1899 IF (DETP <= ZERO) cnvflg = .TRUE.
1900
1901 continue
1902
1903
1904
1905
1906 if (ep_wfn) then
1907 IF ((qw00 .eq. 0.0 .and. qi00 .eq. 0.0)) RETURN
1908 if (ii .eq. 0) then
1909 ii = 1
1910 if (clp .gt. 0.0 .and. clp .lt. 1.0) then
1911 hst(kd) = hstkd
1912 qst(kd) = qstkd
1913 ltl(kd) = ltlkd
1914 q0u(kd) = q0ukd
1915 q0d(kd) = q0dkd
1916 dlb(kd) = dlbkd
1917 qrb(kd) = qrbkd
1918 endif
1919 do l=kd,kb1
1920 FCO(L) = FCO(L) - q0u(l) - q0d(l)
1921 RNN(L) = RNN(L) - q0u(l)*zet(l) - q0d(l)*zet(l+1)
1922 GMH(L) = GMH(L) - q0u(l)*xi(l) - q0d(l)*zet(l+1)
1923 ETA(L) = ZET(L) - ZET(L+1)
1924 GMS(L) = XI(L) - XI(L+1)
1925 Q0U(L) = 0.0
1926 Q0D(L) = 0.0
1927 ENDDO
1928 qw00 = 0.0
1929 qi00 = 0.0
1930
1931
1932
1933
1934 go to 777
1935 else
1936 cnvflg = .true.
1937 endif
1938 endif
1939
1940
1941
1942
1943
1944 = HST(KD) - LTL(KD)*NU*(QST(KD)-QOS)
1945 ST2 = LTL(KD) * VTF(KD)
1946 TEM5 = (QLS + QIS) * eta(kd1)
1947 ST1 = HALF * (TX1-ETA(KD1)*ST1-ST2*(DET-TEM5))*DLB(KD)
1948
1949
1950
1951
1952 = WFN + ST1
1953 AKM = AKM - min(ST1,ZERO)
1954
1955
1956 (KD) = ST1 / (ETA(KD1)*qrb(kd))
1957
1958
1959
1960
1961 = DETP
1962 HCC = HCCP
1963 AKM = AKM / WFN
1964
1965
1966
1967
1968
1969
1970 IF (WRKFUN) THEN
1971 IF (WFN .GE. 0.0) WFNC = WFN
1972 RETURN
1973 ELSEIF (.NOT. CRTFUN) THEN
1974 ACR = WFNC
1975 ENDIF
1976
1977
1978
1979 = .FALSE.
1980
1981 TEM = max(0.05, MIN(CD*200.0, MAX_NEG_BOUY))
1982 IF (WFN .GT. ACR .AND. (.NOT. cnvflg) &
1983
1984 .and. dpneg .lt. 150.0 .AND. AKM .LE. TEM) THEN
1985
1986
1987 = .TRUE.
1988 ENDIF
1989
1990
1991
1992
1993
1994
1995 IF (.NOT. CALCUP) RETURN
1996
1997
1998 IF (ALMIN2 .NE. 0.0) THEN
1999 IF (ALMIN1 .NE. ALMIN2) ST1 = 1.0 / max(ONE_M10,(ALMIN2-ALMIN1))
2000 IF (ALM .LT. ALMIN2) THEN
2001 CLP = CLP * max(0.0, min(1.0,(0.3 + 0.7*(ALM-ALMIN1)*ST1)))
2002
2003
2004 ENDIF
2005 ENDIF
2006
2007
2008
2009 = CLP * RHC
2010 dlq = 0.0
2011 tem = 1.0 / (1.0 + dlq_fac)
2012 do l=kd,kb1
2013 rnn(l) = rns(l) * tem
2014 dlq(l) = rns(l) * tem * dlq_fac
2015 enddo
2016 DO L=KBL,K
2017 RNN(L) = 0.0
2018 ENDDO
2019
2020
2021
2022
2023
2024 = .FALSE.
2025 IF (DNDRFT) THEN
2026
2027 = 0.0
2028 IF (CLP .GT. 0.0) THEN
2029 DO L=KD,KB1
2030 TRAIN = TRAIN + RNN(L)
2031 ENDDO
2032 ENDIF
2033
2034 PL = (PRL(KD1) + PRL(KD))*HALF
2035 TEM = PRL(K+1)*(1.0-DPD*0.001)
2036 IF (TRAIN .GT. 1.0E-4 .AND. PL .LE. TEM) DDFT = .TRUE.
2037
2038 ENDIF
2039
2040
2041
2042
2043
2044
2045
2046 IF (DDFT) THEN
2047 CALL DDRFT( &
2048 & K, KD &
2049 &, TLA, ALFIND &
2050 &, TOL, QOL, HOL, PRL, QST, HST, GAM, GAF &
2051
2052 , QRB, QRT, BUY, KBL, IDH, ETA, RNN, ETAI &
2053 &, ALM, WFN, TRAIN, DDFT &
2054 &, ETD, HOD, QOD, EVP, DOF, CLDFR, ETZ &
2055 &, GMS, GSD, GHD, lprnt)
2056
2057 ENDIF
2058
2059
2060
2061
2062 IF (.NOT. DDFT) THEN
2063 DO L=KD,K+1
2064 ETD(L) = 0.0
2065 HOD(L) = 0.0
2066 QOD(L) = 0.0
2067 ENDDO
2068 DO L=KD,K
2069 EVP(L) = 0.0
2070 ETZ(L) = 0.0
2071 ENDDO
2072
2073 ENDIF
2074
2075
2076
2077
2078
2079
2080
2081
2082 = 0.0
2083
2084
2085
2086
2087
2088
2089 = 0.0
2090
2091
2092
2093 = 1.0 - tem1
2094 TEM = DET * QIL(KD)
2095
2096
2097 st1 = (HCC+ALHF*TEM-ETA(KD)*HST(KD)) / (1.0+gam(KD))
2098 DS = ETA(KD1) * (HOS- HOL(KD)) - ALHL*(QOS - QOL(KD))
2099 DH = ETA(KD1) * (HOS- HOL(KD))
2100
2101
2102 GMS(KD) = (DS + st1 - tem1*det*alhl-tem*alhf) * PRI(KD)
2103 GMH(KD) = PRI(KD) * (HCC-ETA(KD)*HOS + DH)
2104
2105
2106
2107
2108
2109
2110
2111 (KD) = (tem2*(DET-TEM) + ETA(KD1)*(QLS-CLL(KD)) &
2112 & + (1.0-QIL(KD))*dlq(kd) - ETA(KD)*QLS ) * PRI(KD)
2113
2114 QIL(KD) = (tem2*TEM + ETA(KD1)*(QIS-CIL(KD)) &
2115 & + QIL(KD)*dlq(kd) - ETA(KD)*QIS ) * PRI(KD)
2116
2117 (KD) = 0.0
2118 GSD(KD) = 0.0
2119
2120 DO L=KD1,K
2121 ST1 = ONE - ALFINT(L,1)
2122 ST2 = ONE - ALFINT(L,2)
2123 ST3 = ONE - ALFINT(L,3)
2124 ST4 = ONE - ALFINT(L,4)
2125 ST5 = ONE - ALFIND(L)
2126 HB = ALFINT(L,1)*HOL(L-1) + ST1*HOL(L)
2127 QB = ALFINT(L,2)*QOL(L-1) + ST2*QOL(L)
2128
2129 TEM = ALFINT(L,4)*CIL(L-1) + ST4*CIL(L)
2130 TEM2 = ALFINT(L,3)*CLL(L-1) + ST3*CLL(L)
2131
2132 TEM1 = ETA(L) * (TEM - CIL(L))
2133 TEM3 = ETA(L) * (TEM2 - CLL(L))
2134
2135 HBD = ALFIND(L)*HOL(L-1) + ST5*HOL(L)
2136 QBD = ALFIND(L)*QOL(L-1) + ST5*QOL(L)
2137
2138 TEM5 = ETD(L) * (HOD(L) - HBD)
2139 TEM6 = ETD(L) * (QOD(L) - QBD)
2140
2141 = ETA(L) * (HB - HOL(L)) + TEM5
2142 DS = DH - ALHL * (ETA(L) * (QB - QOL(L)) + TEM6)
2143
2144 GMH(L) = DH * PRI(L)
2145 GMS(L) = DS * PRI(L)
2146
2147
2148
2149
2150
2151
2152 (L) = TEM5 * PRI(L)
2153 GSD(L) = (TEM5 - ALHL * TEM6) * PRI(L)
2154
2155 (L) = (TEM3 + (1.0-QIL(L))*dlq(l)) * PRI(L)
2156 QIL(L) = (TEM1 + QIL(L)*dlq(l)) * PRI(L)
2157
2158 TEM1 = ETA(L) * (CIL(L-1) - TEM)
2159 TEM3 = ETA(L) * (CLL(L-1) - TEM2)
2160
2161 DH = ETA(L) * (HOL(L-1) - HB) - TEM5
2162 DS = DH - ALHL * ETA(L) * (QOL(L-1) - QB) &
2163 & + ALHL * (TEM6 - EVP(L-1))
2164
2165 GMH(L-1) = GMH(L-1) + DH * PRI(L-1)
2166 GMS(L-1) = GMS(L-1) + DS * PRI(L-1)
2167
2168
2169
2170
2171
2172 (L-1) = GHD(L-1) - TEM5 * PRI(L-1)
2173 GSD(L-1) = GSD(L-1) - (TEM5-ALHL*(TEM6-EVP(L-1))) * PRI(L-1)
2174
2175 QIL(L-1) = QIL(L-1) + TEM1 * PRI(L-1)
2176 QLL(L-1) = QLL(L-1) + TEM3 * PRI(L-1)
2177
2178
2179
2180
2181
2182
2183 = avh + gmh(l-1)*(prs(l)-prs(l-1))
2184
2185 ENDDO
2186
2187 = HOL(K)
2188 QBD = QOL(K)
2189 TEM5 = ETD(K+1) * (HOD(K+1) - HBD)
2190 TEM6 = ETD(K+1) * (QOD(K+1) - QBD)
2191 DH = - TEM5
2192 DS = DH + ALHL * TEM6
2193 TEM1 = DH * PRI(K)
2194 TEM2 = (DS - ALHL * EVP(K)) * PRI(K)
2195 GMH(K) = GMH(K) + TEM1
2196 GMS(K) = GMS(K) + TEM2
2197 GHD(K) = GHD(K) + TEM1
2198 GSD(K) = GSD(K) + TEM2
2199
2200
2201
2202
2203 = avh + gmh(K)*(prs(KP1)-prs(K))
2204
2205 = - GRAVFAC * pris
2206 TX1 = DH * tem4
2207 TX2 = DS * tem4
2208
2209 DO L=KBL,K
2210 GMH(L) = GMH(L) + TX1
2211 GMS(L) = GMS(L) + TX2
2212 GHD(L) = GHD(L) + TX1
2213 GSD(L) = GSD(L) + TX2
2214
2215 = avh + tx1*(prs(l+1)-prs(l))
2216 ENDDO
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231 DO L=KD,K
2232
2233 TEM1 = GMH(L)
2234 TEM2 = GMS(L)
2235 HOL(L) = HOL(L) + TEM1*TESTMB
2236 QOL(L) = QOL(L) + (TEM1-TEM2) * (TESTMB/ALHL)
2237 HST(L) = HST(L) + TEM2*(ONE+GAM(L))*TESTMB
2238 QST(L) = QST(L) + TEM2*GAM(L)*(TESTMB/ALHL)
2239 CLL(L) = CLL(L) + QLL(L) * TESTMB
2240 CIL(L) = CIL(L) + QIL(L) * TESTMB
2241 ENDDO
2242
2243 if (alm .gt. 0.0) then
2244 HOS = HOS + GMH(KD) * TESTMB
2245 QOS = QOS + (GMH(KD)-GMS(KD)) * (TESTMB/ALHL)
2246 QLS = QLS + QLL(KD) * TESTMB
2247 QIS = QIS + QIL(KD) * TESTMB
2248 else
2249 st2 = 1.0 - st1s
2250 HOS = HOS + (st1s*GMH(KD)+st2*GMH(KD1)) * TESTMB
2251 QOS = QOS + (st1s * (GMH(KD)-GMS(KD)) &
2252 & + st2 * (GMH(KD1)-GMS(KD1))) * (TESTMB/ALHL)
2253 HST(kd) = HST(kd) + (st1s*GMS(kd)*(ONE+GAM(kd)) &
2254 & + st2*gms(kd1)*(ONE+GAM(kd1))) * TESTMB
2255 QST(kd) = QST(kd) + (st1s*GMS(kd)*GAM(kd) &
2256 & + st2*gms(kd1)*gam(kd1)) * (TESTMB/ALHL)
2257
2258 QLS = QLS + (st1s*QLL(KD)+st2*QLL(KD1)) * TESTMB
2259 QIS = QIS + (st1s*QIL(KD)+st2*QIL(KD1)) * TESTMB
2260 endif
2261
2262
2263 = PRL(Kmaxp1) - PRL(Kmax)
2264 HBL = HOL(Kmax) * TEM
2265 QBL = QOL(Kmax) * TEM
2266 QLB = CLL(Kmax) * TEM
2267 QIB = CIL(Kmax) * TEM
2268 DO L=KmaxM1,KBL,-1
2269 TEM = PRL(L+1) - PRL(L)
2270 HBL = HBL + HOL(L) * TEM
2271 QBL = QBL + QOL(L) * TEM
2272 QLB = QLB + CLL(L) * TEM
2273 QIB = QIB + CIL(L) * TEM
2274 ENDDO
2275 HBL = HBL * PRISM
2276 QBL = QBL * PRISM
2277 QLB = QLB * PRISM
2278 QIB = QIB * PRISM
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291 = ZERO
2292 TX1 = ZERO
2293 QTL = QST(KB1) - GAF(KB1)*HST(KB1)
2294 QTV = QBL
2295 HCC = HBL
2296 TX2 = HCC
2297 TX4 = (ALHF*0.5)*MAX(ZERO,MIN(ONE,(TCR-TCL-TOL(KB1))*TCRF))
2298
2299 = qbl
2300 tx1 = qib + qlb
2301
2302
2303 DO L=KB1,KD1,-1
2304 DEL_ETA = ETA(L) - ETA(L+1)
2305 HCCP = HCC + DEL_ETA*HOL(L)
2306
2307 = QST(L-1) - GAF(L-1)*HST(L-1)
2308 QTVP = 0.5 * ((QTLP+QTL)*ETA(L) + (GAF(L)+GAF(L-1))*HCCP)
2309
2310 DETP = (BKC(L)*TX1 - (QTVP-QTV) &
2311 & + DEL_ETA*(QOL(L)+CLL(L)+CIL(L)) &
2312 & + ETA(L)*Q0U(L) + ETA(L+1)*Q0D(L)) * AKC(L)
2313 IF (DETP .LE. ZERO) cnvflg = .TRUE.
2314
2315 ST1 = HST(L) - LTL(L)*NU*(QST(L)-QOL(L))
2316
2317 TEM2 = (ALHF*0.5)*MAX(ZERO,MIN(ONE,(TCR-TCL-TOL(L-1))*TCRF))
2318 TEM1 = HCCP + DETP * (TEM2+TX4)
2319
2320 ST2 = LTL(L) * VTF(L)
2321 TEM5 = CLL(L) + CIL(L)
2322 AKM = AKM + &
2323 & ( (TX2 -ETA(L+1)*ST1-ST2*(TX1-TEM5*eta(l+1))) * DLB(L) &
2324 & + (TEM1 -ETA(L )*ST1-ST2*(DETP-TEM5*eta(l))) * DLT(L) )
2325
2326 = HCCP
2327 TX1 = DETP
2328 TX2 = TEM1
2329 QTL = QTLP
2330 QTV = QTVP
2331 TX4 = TEM2
2332 ENDDO
2333
2334 if (cnvflg) return
2335
2336
2337
2338
2339
2340
2341 = HST(KD) - LTL(KD)*NU*(QST(KD)-QOS)
2342 ST2 = LTL(KD) * VTF(KD)
2343 TEM5 = (QLS + QIS) * eta(kd1)
2344 AKM = AKM + HALF * (TX2-ETA(KD1)*ST1-ST2*(TX1-TEM5)) * DLB(KD)
2345
2346 = (AKM - WFN) * (ONE/TESTMB)
2347
2348
2349
2350
2351
2352
2353 = rel_fac
2354
2355 = - (WFN-ACR) / AKM
2356
2357
2358
2359
2360
2361
2362
2363 = AMB * CLP * tem2
2364
2365
2366
2367
2368
2369 = (PRL(KMAXP1)-PRL(KBL))*(FRACBL*GRAVCON)
2370 AMB = MAX(MIN(AMB, AMBMAX),ZERO)
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380 = 0.0
2381 avq = 0.0
2382 avr = dof
2383
2384
2385 = DSFC + AMB * ETD(K) * (1.0/DT)
2386
2387
2388 DO L=K,KD,-1
2389 PCU(L) = PCU(L) + AMB*RNN(L)
2390 = avr + rnn(l)
2391
2392 ENDDO
2393
2394
2395
2396 = AMB * (ONE/CP)
2397 TX2 = AMB * (ONE/ALHL)
2398 DO L=KD,K
2399 ST1 = GMS(L)*TX1
2400 TOI(L) = TOI(L) + ST1
2401 TCU(L) = TCU(L) + ST1
2402 TCD(L) = TCD(L) + GSD(L) * TX1
2403
2404 = st1 - (alhl/cp) * (QIL(L) + QLL(L)) * AMB
2405
2406 avt = avt + st1 * (prs(l+1)-prs(l))
2407
2408 FLX(L) = FLX(L) + ETA(L)*AMB
2409 FLXD(L) = FLXD(L) + ETD(L)*AMB
2410
2411 (L) = QII(L) + QIL(L) * AMB
2412 TEM = 0.0
2413
2414 QLI(L) = QLI(L) + QLL(L) * AMB + TEM
2415
2416 ST1 = (GMH(L)-GMS(L)) * TX2
2417
2418 QOI(L) = QOI(L) + ST1
2419 QCU(L) = QCU(L) + ST1
2420 QCD(L) = QCD(L) + (GHD(L)-GSD(L)) * TX2
2421
2422 = avq + (st1+(QLL(L)+QIL(L))*amb) * (prs(l+1)-prs(l))
2423
2424
2425
2426
2427
2428
2429
2430
2431 ENDDO
2432 avr = avr * amb
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447
2448
2449
2450
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464
2465
2466
2467
2468
2469
2470 = 0.0
2471 TX2 = 0.0
2472
2473 IF (REVAP) THEN
2474
2475 = 0.0
2476 do l=kd,kbl
2477 IF (L .lt. IDH .or. (.not. DDFT)) THEN
2478 tem = tem + amb * rnn(l)
2479 endif
2480 enddo
2481 tem = tem + amb * dof
2482 tem = tem * (3600.0/dt)
2483
2484
2485
2486
2487 = sqrt(max(1.0, min(100.0,(4.0E10/max(garea,one)))))
2488
2489
2490
2491 = max(ZERO, min(ONE, rknob*clf(tem)*tem1))
2492
2493
2494
2495
2496
2497
2498
2499
2500 DO L=KD,KBL
2501
2502 IF (L .GE. IDH .AND. DDFT) THEN
2503 TX2 = TX2 + AMB * RNN(L)
2504 CLDFRD = MIN(AMB*CLDFR(L), clfrac)
2505 ELSE
2506 TX1 = TX1 + AMB * RNN(L)
2507 ENDIF
2508 tx4 = zfac * phil(l)
2509 tx4 = (one - tx4 * (one - half*tx4)) * afc
2510
2511 IF (TX1 .GT. 0. .OR. TX2 .GT. 0.0) THEN
2512 TEQ = TOI(L)
2513 QEQ = QOI(L)
2514 PL = 0.5 * (PRL(L+1)+PRL(L))
2515
2516 ST1 = MAX(ZERO, MIN(ONE, (TCR-TEQ)*TCRF))
2517 ST2 = ST1*ELFOCP + (1.0-ST1)*ELOCP
2518
2519 CALL QSATCN ( TEQ,PL,QSTEQ,DQDT)
2520
2521
2522 = 0.5 * (QSTEQ*rhc_ls(l)-QEQ) / (1.+ST2*DQDT)
2523
2524 = QEQ + DELTAQ
2525 TEQ = TEQ - DELTAQ*ST2
2526
2527 = MAX(ZERO, MIN(ONE, (TCR-TEQ)*TCRF))
2528 TEM2 = TEM1*ELFOCP + (1.0-TEM1)*ELOCP
2529
2530 CALL QSATCN ( TEQ,PL,QSTEQ,DQDT)
2531
2532
2533 = (QSTEQ*rhc_ls(l)-QEQ) / (1.+TEM2*DQDT)
2534
2535 = QEQ + DELTAQ
2536 TEQ = TEQ - DELTAQ*TEM2
2537
2538 IF (QEQ .GT. QOI(L)) THEN
2539 POTEVAP = (QEQ-QOI(L))*(PRL(L+1)-PRL(L))*GRAVCON
2540
2541 tem4 = 0.0
2542 if (tx1 .gt. 0.0) &
2543 & TEM4 = POTEVAP * (1. - EXP( tx4*TX1**0.57777778 ) )
2544
2545 = MIN(TX1, TEM4*CLFRAC)
2546
2547
2548
2549
2550
2551
2552 if (tx1 .lt. rainmin*dt) actevap = min(tx1, potevap)
2553
2554 = 0.0
2555 if (tx2 .gt. 0.0) &
2556 & TEM4 = POTEVAP * (1. - EXP( tx4*TX2**0.57777778 ) )
2557
2558 = min(MIN(TX2, TEM4*CLDFRD), potevap-actevap)
2559 if (tx2 .lt. rainmin*dt) tem4 = min(tx2, potevap-actevap)
2560
2561 = TX1 - ACTEVAP
2562 TX2 = TX2 - TEM4
2563 ST1 = (ACTEVAP+TEM4) * PRI(L)
2564 QOI(L) = QOI(L) + ST1
2565 QCU(L) = QCU(L) + ST1
2566
2567
2568 = ST1 * ELOCP
2569 TOI(L) = TOI(L) - ST1
2570 TCU(L) = TCU(L) - ST1
2571 ENDIF
2572 ENDIF
2573 ENDDO
2574
2575 = CUP + TX1 + TX2 + DOF * AMB
2576 ELSE
2577 DO L=KD,K
2578 TX1 = TX1 + AMB * RNN(L)
2579 ENDDO
2580 CUP = CUP + TX1 + DOF * AMB
2581 ENDIF
2582
2583
2584
2585
2586
2587
2588
2589
2590 if (NTRC > 0) then
2591 do l=kd,k-1
2592 if (etz(l) .ne. zero) etzi(l) = one / etz(l)
2593 enddo
2594 DO N=1,NTRC
2595
2596 DO L=KD,K
2597 HOL(L) = ROI(L,N)
2598 ENDDO
2599
2600 = RBL(N)
2601 HOD(KD) = HOL(KD)
2602
2603 DO L=KD1,K
2604 ST1 = ONE - ALFIND(L)
2605 HB = ALFIND(L) * HOL(L-1) + ST1 * HOL(L)
2606 IF (ETZ(L-1) .NE. ZERO) THEN
2607 TEM = ETZI(L-1)
2608 IF (ETD(L) > ETD(L-1)) THEN
2609 HOD(L) = (ETD(L-1)*(HOD(L-1)-HOL(L-1)) &
2610 & + ETD(L) *(HOL(L-1)-HB) + ETZ(L-1)*HB) * TEM
2611 ELSE
2612 HOD(L) = (ETD(L-1)*(HOD(L-1)-HB) + ETZ(L-1)*HB) * TEM
2613 ENDIF
2614 ELSE
2615 HOD(L) = HB
2616 ENDIF
2617 ENDDO
2618
2619 DO L=KB1,KD,-1
2620 HCC = HCC + (ETA(L)-ETA(L+1))*HOL(L)
2621 ENDDO
2622
2623
2624
2625
2626
2627
2628 if (FSCAV_(N) > 0.0) then
2629 DELZKM = ( PHIL(KD) - PHIH(KD1) ) *(onebg*0.001)
2630 FNOSCAV = exp(- FSCAV_(N) * DELZKM)
2631 else
2632 FNOSCAV = 1.0
2633 endif
2634
2635 GMH(KD) = PRI(KD) * (HCC-ETA(KD)*HOL(KD)) * trcfac(kd,n) &
2636 & * FNOSCAV
2637 DO L=KD1,K
2638 if (FSCAV_(N) > 0.0) then
2639 DELZKM = ( PHIL(KD) - PHIH(L+1) ) *(onebg*0.001)
2640 FNOSCAV = exp(- FSCAV_(N) * DELZKM)
2641 endif
2642 ST1 = ONE - ALFINT(L,N+4)
2643 ST2 = ONE - ALFIND(L)
2644 HB = ALFINT(L,N+4) * HOL(L-1) + ST1 * HOL(L)
2645 HBD = ALFIND(L) * HOL(L-1) + ST2 * HOL(L)
2646 TEM5 = ETD(L) * (HOD(L) - HBD)
2647 DH = ETA(L) * (HB - HOL(L)) * FNOSCAV + TEM5
2648 GMH(L ) = DH * PRI(L) * trcfac(l,n)
2649 DH = ETA(L) * (HOL(L-1) - HB) * FNOSCAV - TEM5
2650 GMH(L-1) = GMH(L-1) + DH * PRI(L-1) * trcfac(l,n)
2651 ENDDO
2652
2653 DO L=KD,K
2654 ST1 = GMH(L)*AMB
2655 ROI(L,N) = HOL(L) + ST1
2656 RCU(L,N) = RCU(L,N) + ST1
2657 ENDDO
2658 ENDDO
2659 endif
2660
2661
2662
2663
2664 RETURN
2665 END
2666
2667 SUBROUTINE DDRFT( &
2668 & K, KD &
2669 &, TLA, ALFIND &
2670 &, TOL, QOL, HOL, PRL, QST, HST, GAM, GAF &
2671
2672 , QRB, QRT, BUY, KBL, IDH, ETA, RNN, ETAI &
2673 &, ALM, WFN, TRAIN, DDFT &
2674 &, ETD, HOD, QOD, EVP, DOF, CLDFRD, WCB &
2675 &, GMS, GSD, GHD,lprnt)
2676
2677
2678
2679
2680
2681
2682
2683
2684
2685
2686
2687
2688
2689
2690
2691
2692
2693
2694
2695
2696
2697
2698 USE MACHINE , ONLY : kind_phys
2699 use module_ras
2700 IMPLICIT NONE
2701
2702
2703
2704 INTEGER K, KD
2705 real(kind=kind_phys) ALFIND(K)
2706 INTEGER KBL, KB1
2707
2708 LOGICAL SKPDD, SKPUP
2709
2710 real(kind=kind_phys) HOL(KD:K), QOL(KD:K), GAF(KD:K+1) &
2711 &, HST(KD:K), QST(KD:K), TOL(KD:K) &
2712 &, BUY(KD:K+1), QRB(KD:K), QRT(KD:K) &
2713 &, GAM(KD:K+1), RNN(KD:K), RNS(KD:K) &
2714 &, ETA(KD:K+1), PRL(KD:K+1), ETAI(KD:K)
2715
2716
2717
2718 real(kind=kind_phys) TRAIN, WFN, ALM
2719
2720
2721
2722 real(kind=kind_phys) GMS(KD:K+1)
2723 real(kind=kind_phys) TX1, TX2, TX3, TX4 &
2724 &, TX5, TX6, TX7, TX8, TX9
2725 LOGICAL cnvflg
2726
2727 real(kind=kind_phys) TL, PL, QL, QS, DQS, ST1, HB, QB, TB &
2728 &, QQQ, PICON, PIINV, DEL_ETA &
2729 &, TEM, TEM1, TEM2, TEM3, TEM4, ST2 &
2730 &, ERRMIN, ERRMI2, ERRH, ERRW, ERRE, TEM5 &
2731 &, TEM6, HBD, QBD
2732 INTEGER I, L, N, IX, KD1, II &
2733 &, KP1, IT, KM1, KTEM, KK, KK1, LM1, LL, LP1 &
2734 &, IP1, JJ, ntla
2735
2736
2737 integer, parameter :: NUMTLA=2
2738
2739 parameter (ERRMIN=0.0001, ERRMI2=0.1*ERRMIN)
2740
2741
2742 real(kind=kind_phys) TLA, STLA, CTL2, CTL3
2743 real(kind=kind_phys) GMF, PI, ONPG, CTLA, VTRM, VTPEXP &
2744 &, RPART, QRMIN, AA1, BB1, CC1, DD1 &
2745 &, WC2MIN, WCMIN, WCBASE, F2, F3, F5, GMF1, GMF5 &
2746 &, QRAF, QRBF, del_tla
2747
2748
2749 parameter (ONPG=1.0+0.5, GMF=1.0/ONPG, RPART=0.0)
2750
2751
2752
2753
2754 PARAMETER (AA1=1.0, BB1=1.0, CC1=1.0, DD1=1.0, F3=CC1, F5=1.0)
2755 parameter (QRMIN=1.0E-6, WC2MIN=0.01, GMF1=GMF/AA1, GMF5=GMF/F5)
2756
2757
2758
2759 PARAMETER (PI=3.1415926535897931, PIINV=1.0/PI)
2760 INTEGER ITR, ITRMU, ITRMD, KTPD, ITRMIN, ITRMND
2761
2762 PARAMETER (ITRMU=25, ITRMD=25, ITRMIN=12, ITRMND=12)
2763
2764
2765
2766 real(kind=kind_phys) QRP(KD:K+1), WVL(KD:K+1), AL2
2767 real(kind=kind_phys) WVLO(KD:K+1)
2768
2769 real(kind=kind_phys) RNF(KD:K), ETD(KD:K+1), WCB(KD:K) &
2770 &, HOD(KD:K+1), QOD(KD:K+1), EVP(KD:K) &
2771 &, ROR(KD:K+1), STLT(KD:K) &
2772 &, GHD(KD:K), GSD(KD:K), CLDFRD(KD:K) &
2773 &, RNT, RNB &
2774 &, ERRQ, RNTP
2775 INTEGER IDW, IDH, IDN(K), idnm
2776 real(kind=kind_phys) ELM(K)
2777
2778 real(kind=kind_phys) EDZ, DDZ, CE, QHS, FAC, FACG, ASIN, &
2779 & RSUM1, RSUM2, RSUM3, CEE
2780 LOGICAL DDFT, UPDRET, DDLGK
2781
2782 real(kind=kind_phys) AA(KD:K,KD:K+1), QW(KD:K,KD:K) &
2783 &, BUD(KD:K), VT(2), VRW(2), TRW(2) &
2784 &, GQW(KD:K) &
2785 &, QA(3), WA(3), DOF, DOFW &
2786 &, QRPI(KD:K), QRPS(KD:K)
2787
2788
2789
2790
2791 real(kind=kind_phys) QRPF, VTPF
2792 logical lprnt
2793
2794
2795
2796 = KD + 1
2797 KP1 = K + 1
2798 KM1 = K - 1
2799 KB1 = KBL - 1
2800
2801
2802 = -0.3636
2803
2804 = PI * ONEBG * 0.5
2805
2806
2807
2808 = 0.0
2809 RNTP = 0.0
2810 DOF = 0.0
2811 ERRQ = 10.0
2812 RNB = 0.0
2813 RNT = 0.0
2814 TX2 = PRL(KBL)
2815
2816 = (PRL(KD) + PRL(KD1)) * 0.5
2817 ROR(KD) = CMPOR*TX1 / (TOL(KD)*(1.0+NU*QOL(KD)))
2818
2819 (KD) = VTP * VTPF(ROR(KD))
2820
2821 (KD) = QRMIN
2822
2823 = TOL(K) * (1.0 + NU * QOL(K))
2824 ROR(K+1) = 0.5 * CMPOR * (PRL(K+1)+PRL(K)) / TEM
2825 GMS(K+1) = VTP * VTPF(ROR(K+1))
2826 QRP(K+1) = QRMIN
2827
2828 = kbl
2829 DO L=KD1,K
2830 TEM = 0.5 * (TOL(L)+TOL(L-1)) &
2831 & * (1.0 + (0.5*NU) * (QOL(L)+QOL(L-1)))
2832 ROR(L) = CMPOR * PRL(L) / TEM
2833
2834 (L) = VTP * VTPF(ROR(L))
2835 QRP(L) = QRMIN
2836 if (buy(l) .le. 0.0 .and. kk .eq. KBL) then
2837 kk = l
2838 endif
2839 ENDDO
2840 if (kk .ne. kbl) then
2841 do l=kk,kbl
2842 buy(l) = 0.9 * buy(l-1)
2843 enddo
2844 endif
2845
2846 do l=kd,k
2847 qrpi(l) = buy(l)
2848 enddo
2849 do l=kd1,kb1
2850 buy(l) = 0.25 * (qrpi(l-1)+qrpi(l)+qrpi(l)+qrpi(l+1))
2851 enddo
2852
2853
2854
2855 = 1000.0 + tx1 - prl(k+1)
2856 CALL ANGRAD(TX1, ALM, AL2, TLA, TX2, WFN, TX3)
2857
2858
2859
2860 = 2.0*BB1*ONEBG/(PI*0.2)
2861 WCMIN = SQRT(WC2MIN)
2862 WCBASE = WCMIN
2863
2864
2865
2866 = TLA * 0.3
2867 TLA = TLA - DEL_TLA
2868
2869 DO L=KD,K
2870 RNF(L) = 0.0
2871 RNS(L) = 0.0
2872 WVL(L) = 0.0
2873 STLT(L) = 0.0
2874 GQW(L) = 0.0
2875 QRP(L) = QRMIN
2876 DO N=KD,K
2877 QW(N,L) = 0.0
2878 ENDDO
2879 ENDDO
2880
2881
2882
2883 = KBL
2884 QW(KD,KD) = -QRB(KD) * GMF1
2885 GHD(KD) = ETA(KD) * ETA(KD)
2886 GQW(KD) = QW(KD,KD) * GHD(KD)
2887 GSD(KD) = ETAI(KD) * ETAI(KD)
2888
2889 (KK) = - QRB(KK-1) * (GMF1+GMF1)
2890
2891 (KK) = WCBASE * WCBASE
2892
2893 TX1 = WCB(KK)
2894 GSD(KK) = 1.0
2895 GHD(KK) = 1.0
2896
2897 = GMF1 + GMF1
2898 DO L=KB1,KD1,-1
2899 GHD(L) = ETA(L) * ETA(L)
2900 GSD(L) = ETAI(L) * ETAI(L)
2901 GQW(L) = - GHD(L) * (QRB(L-1)+QRT(L)) * TEM
2902 QW(L,L) = - QRT(L) * TEM
2903
2904 = 0.5 * (eta(l) + eta(l+1))
2905 TX1 = TX1 + BUY(L) * TEM * (qrb(l)+qrt(l)) * st1 * st1
2906 WCB(L) = TX1 * GSD(L)
2907 ENDDO
2908
2909 = (QRB(KD) + QRT(KD1) + QRT(KD1)) * GMF1
2910 GQW(KD1) = - GHD(KD1) * TEM1
2911 QW(KD1,KD1) = - QRT(KD1) * TEM
2912 st1 = 0.5 * (eta(kd) + eta(kd1))
2913 WCB(KD) = (TX1 + BUY(KD)*TEM*qrb(kd)*st1*st1) * GSD(KD)
2914
2915 DO L=KD1,KBL
2916 DO N=KD,L-1
2917 QW(N,L) = GQW(L) * GSD(N)
2918 ENDDO
2919 ENDDO
2920 QW(KBL,KBL) = 0.0
2921
2922 do ntla=1,numtla
2923
2924
2925 if (errq .lt. 0.1 .or. tla .gt. 45.0) cycle
2926
2927 = tla + del_tla
2928 STLA = SIN(TLA*PI/180.0)
2929 CTL2 = 1.0 - STLA * STLA
2930
2931
2932
2933
2934
2935 = F2 * STLA * AL2
2936 CTL2 = DD1 * CTL2
2937 CTL3 = 0.1364 * CTL2
2938
2939 DO L=KD,K
2940 RNF(L) = 0.0
2941 WVL(L) = 0.0
2942 STLT(L) = 0.0
2943 QRP(L) = QRMIN
2944 ENDDO
2945 WVL(KBL) = WCBASE
2946 STLT(KBL) = 1.0 / WCBASE
2947
2948 DO L=KD,K+1
2949 DO N=KD,K
2950 AA(N,L) = 0.0
2951 ENDDO
2952 ENDDO
2953
2954 = .FALSE.
2955
2956 DO ITR=1,ITRMU
2957 IF (.NOT. SKPUP) THEN
2958 wvlo = wvl
2959
2960
2961
2962 = 0.0
2963 QRPI(KBL) = 1.0 / QRP(KBL)
2964 DO L=KB1,KD,-1
2965 TX1 = TX1 + QRP(L+1) * GQW(L+1)
2966 ST1 = WCB(L) + QW(L,L) * QRP(L) &
2967 & + TX1 * GSD(L)
2968 if (st1 .gt. wc2min) then
2969
2970 (L) = 0.5 * (SQRT(ST1) + WVL(L))
2971
2972 else
2973
2974
2975
2976
2977
2978 (l) = 0.5*(wvl(l) + wvl(l+1))
2979 qrp(l) = 0.5*((wvl(l)*wvl(l)-wcb(l)-tx1*gsd(l))/qw(l,l) &
2980 & + qrp(l))
2981
2982 endif
2983
2984
2985 (l) = max(wvl(l), wcbase)
2986 STLT(L) = 1.0 / WVL(L)
2987 QRPI(L) = 1.0 / QRP(L)
2988 ENDDO
2989
2990
2991
2992
2993
2994
2995
2996
2997
2998
2999
3000
3001 (1) = GMS(KD) * QRPF(QRP(KD))
3002 TRW(1) = ETA(KD) * QRP(KD) * STLT(KD)
3003 TX6 = TRW(1) * VT(1)
3004 VRW(1) = F3*WVL(KD) - CTL2*VT(1)
3005 BUD(KD) = STLA * TX6 * QRB(KD) * 0.5
3006 RNF(KD) = BUD(KD)
3007 DOF = 1.1364 * BUD(KD) * QRPI(KD)
3008 DOFW = -BUD(KD) * STLT(KD)
3009
3010 = TRW(1) * VRW(1)
3011 TX2 = 0.0
3012 TX4 = 0.0
3013 RNB = RNT
3014 TX1 = 0.5
3015 TX8 = 0.0
3016
3017 IF (RNT .GE. 0.0) THEN
3018 TX3 = (RNT-CTL3*TX6) * QRPI(KD)
3019 TX5 = CTL2 * TX6 * STLT(KD)
3020 ELSE
3021 TX3 = 0.0
3022 TX5 = 0.0
3023 RNT = 0.0
3024 RNB = 0.0
3025 ENDIF
3026
3027 DO L=KD1,KB1
3028 KTEM = MAX(L-2, KD)
3029 LL = L - 1
3030
3031
3032 (2) = GMS(L) * QRPF(QRP(L))
3033 TRW(2) = ETA(L) * QRP(L) * STLT(L)
3034 VRW(2) = F3*WVL(L) - CTL2*VT(2)
3035 QQQ = STLA * TRW(2) * VT(2)
3036 ST1 = TX1 * QRB(LL)
3037 BUD(L) = QQQ * (ST1 + QRT(L))
3038
3039 (2) = DOF
3040 WA(2) = DOFW
3041 DOF = 1.1364 * BUD(L) * QRPI(L)
3042 DOFW = -BUD(L) * STLT(L)
3043
3044 (LL) = RNF(LL) + QQQ * ST1
3045 RNF(L) = QQQ * QRT(L)
3046
3047 = VRW(1) + VRW(2)
3048 TEM4 = TRW(1) + TRW(2)
3049
3050 = .25 * TEM3 * TEM4
3051 TEM4 = TEM4 * CTL3
3052
3053
3054
3055
3056 = .25*(TRW(1)*TEM3 - TEM4*VT(1))*QRPI(LL)
3057 ST1 = .25*(TRW(1)*(CTL2*VT(1)-VRW(2)) &
3058 & * STLT(LL) + F3*TRW(2))
3059
3060 = .25*(TRW(2)*TEM3 - TEM4*VT(2))*QRPI(L)
3061 ST2 = .25*(TRW(2)*(CTL2*VT(2)-VRW(1)) &
3062 & * STLT(L) + F3*TRW(1))
3063
3064
3065
3066 (1) = TX2
3067 QA(2) = QA(2) + TX3 - TEM1
3068 QA(3) = -TEM2
3069
3070 (1) = TX4
3071 WA(2) = WA(2) + TX5 - ST1
3072 WA(3) = -ST2
3073
3074 = TEM1
3075 TX3 = TEM2
3076 TX4 = ST1
3077 TX5 = ST2
3078
3079 (1) = VT(2)
3080 TRW(1) = TRW(2)
3081 VRW(1) = VRW(2)
3082
3083 IF (WVL(KTEM) .EQ. WCMIN) WA(1) = 0.0
3084 IF (WVL(LL) .EQ. WCMIN) WA(2) = 0.0
3085 IF (WVL(L) .EQ. WCMIN) WA(3) = 0.0
3086 DO N=KTEM,KBL
3087 AA(LL,N) = (WA(1)*QW(KTEM,N) * STLT(KTEM) &
3088 & + WA(2)*QW(LL,N) * STLT(LL) &
3089 & + WA(3)*QW(L,N) * STLT(L) ) * 0.5
3090 ENDDO
3091 AA(LL,KTEM) = AA(LL,KTEM) + QA(1)
3092 AA(LL,LL) = AA(LL,LL) + QA(2)
3093 AA(LL,L) = AA(LL,L) + QA(3)
3094 BUD(LL) = (TX8 + RNN(LL)) * 0.5 &
3095 & - RNB + TX6 - BUD(LL)
3096 AA(LL,KBL+1) = BUD(LL)
3097 RNB = TX6
3098 TX1 = 1.0
3099 TX8 = RNN(LL)
3100 ENDDO
3101 L = KBL
3102 LL = L - 1
3103
3104 (2) = GMS(L) * QRPF(QRP(L))
3105 TRW(2) = ETA(L) * QRP(L) * STLT(L)
3106 VRW(2) = F3*WVL(L) - CTL2*VT(2)
3107 ST1 = STLA * TRW(2) * VT(2) * QRB(LL)
3108 BUD(L) = ST1
3109
3110 QA(2) = DOF
3111 WA(2) = DOFW
3112 DOF = 1.1364 * BUD(L) * QRPI(L)
3113 DOFW = -BUD(L) * STLT(L)
3114
3115 (LL) = RNF(LL) + ST1
3116
3117 = VRW(1) + VRW(2)
3118 TEM4 = TRW(1) + TRW(2)
3119
3120 = .25 * TEM3 * TEM4
3121 TEM4 = TEM4 * CTL3
3122
3123
3124
3125 = .25*(TRW(1)*TEM3 - TEM4*VT(1))*QRPI(LL)
3126 ST1 = .25*(TRW(1)*(CTL2*VT(1)-VRW(2)) &
3127 & * STLT(LL) + F3*TRW(2))
3128
3129 = .25*(TRW(2)*TEM3 - TEM4*VT(2))*QRPI(L)
3130 ST2 = .25*(TRW(2)*(CTL2*VT(2)-VRW(1)) &
3131 & * STLT(L) + F3*TRW(1))
3132
3133
3134
3135 (1) = TX2
3136 QA(2) = QA(2) + TX3 - TEM1
3137 QA(3) = -TEM2
3138
3139 (1) = TX4
3140 WA(2) = WA(2) + TX5 - ST1
3141 WA(3) = -ST2
3142
3143 = TEM1
3144 TX3 = TEM2
3145 TX4 = ST1
3146 TX5 = ST2
3147
3148 = MAX(L-2, KD)
3149
3150 IF (WVL(IDW) .EQ. WCMIN) WA(1) = 0.0
3151 IF (WVL(LL) .EQ. WCMIN) WA(2) = 0.0
3152 IF (WVL(L) .EQ. WCMIN) WA(3) = 0.0
3153
3154 = IDW
3155 DO N=KK,L
3156 AA(LL,N) = (WA(1)*QW(KK,N) * STLT(KK) &
3157 & + WA(2)*QW(LL,N) * STLT(LL) &
3158 & + WA(3)*QW(L,N) * STLT(L) ) * 0.5
3159
3160 ENDDO
3161
3162 (LL,IDW) = AA(LL,IDW) + QA(1)
3163 AA(LL,LL) = AA(LL,LL) + QA(2)
3164 AA(LL,L) = AA(LL,L) + QA(3)
3165 BUD(LL) = (TX8+RNN(LL)) * 0.5 - RNB + TX6 - BUD(LL)
3166
3167 (LL,L+1) = BUD(LL)
3168
3169 = TRW(2) * VRW(2)
3170
3171
3172
3173 IF (RNB .LT. 0.0) THEN
3174 KK = KBL
3175 TEM = VT(2) * TRW(2)
3176 QA(2) = (RNB - CTL3*TEM) * QRPI(KK)
3177 WA(2) = CTL2 * TEM * STLT(KK)
3178 ELSE
3179 RNB = 0.0
3180 QA(2) = 0.0
3181 WA(2) = 0.0
3182 ENDIF
3183
3184 (1) = TX2
3185 QA(2) = DOF + TX3 - QA(2)
3186 QA(3) = 0.0
3187
3188 (1) = TX4
3189 WA(2) = DOFW + TX5 - WA(2)
3190 WA(3) = 0.0
3191
3192 = KBL
3193 IF (WVL(KK-1) .EQ. WCMIN) WA(1) = 0.0
3194 IF (WVL(KK) .EQ. WCMIN) WA(2) = 0.0
3195
3196 DO II=1,2
3197 N = KK + II - 2
3198 AA(KK,N) = (WA(1)*QW(KK-1,N) * STLT(KK-1) &
3199 & + WA(2)*QW(KK,N) * STLT(KK)) * 0.5
3200 ENDDO
3201 FAC = 0.5
3202 LL = KBL
3203 L = LL + 1
3204 LM1 = LL - 1
3205 AA(LL,LM1) = AA(LL,LM1) + QA(1)
3206 AA(LL,LL) = AA(LL,LL) + QA(2)
3207 BUD(LL) = 0.5*RNN(LM1) - TX6 + RNB - BUD(LL)
3208 AA(LL,LL+1) = BUD(LL)
3209
3210
3211
3212 DO L=KD1,KBL
3213 LM1 = L - 1
3214 cnvflg = ABS(AA(LM1,LM1)) .LT. ABS(AA(L,LM1))
3215 DO N=LM1,KBL+1
3216 IF (cnvflg) THEN
3217 TX1 = AA(LM1,N)
3218 AA(LM1,N) = AA(L,N)
3219 AA(L,N) = TX1
3220 ENDIF
3221 ENDDO
3222 TX1 = AA(L,LM1) / AA(LM1,LM1)
3223 DO N=L,KBL+1
3224 AA(L,N) = AA(L,N) - TX1 * AA(LM1,N)
3225 ENDDO
3226 ENDDO
3227
3228
3229
3230 = KBL
3231 KK1 = KK + 1
3232 AA(KK,KK1) = AA(KK,KK1) / AA(KK,KK)
3233 = ABS(AA(KK,KK1)) * QRPI(KK)
3234
3235
3236
3237 = KBL + 1
3238 DO L=KB1,KD,-1
3239 LP1 = L + 1
3240 TX1 = 0.0
3241 DO N=LP1,KBL
3242 TX1 = TX1 + AA(L,N) * AA(N,KK)
3243 ENDDO
3244 AA(L,KK) = (AA(L,KK) - TX1) / AA(L,L)
3245 = MAX(TX2, ABS(AA(L,KK))*QRPI(L))
3246
3247
3248
3249
3250 ENDDO
3251
3252
3253 if (tx2 .gt. 1.0 .and. abs(errq-tx2) .gt. 0.1) then
3254 tem = 0.5
3255
3256
3257 else
3258 tem = 1.0
3259 endif
3260
3261 DO L=KD,KBL
3262
3263 (L) = MAX(QRP(L)+AA(L,KBL+1)*tem, QRMIN)
3264 ENDDO
3265
3266
3267
3268 IF (ITR .LT. ITRMIN) THEN
3269 TEM = ABS(ERRQ-TX2)
3270 IF (TEM .GE. ERRMI2 .AND. TX2 .GE. ERRMIN) THEN
3271 ERRQ = TX2
3272 ELSE
3273 SKPUP = .TRUE.
3274 = 0.0
3275
3276
3277 ENDIF
3278 ELSE
3279 TEM = ERRQ - TX2
3280
3281 IF (TEM .LT. ZERO .AND. ERRQ .GT. 0.5) THEN
3282
3283
3284
3285 = .TRUE.
3286 = 10.0
3287
3288 ELSEIF (TX2.LT.ERRMIN) THEN
3289 SKPUP = .TRUE.
3290 = 0.0
3291
3292 elseif (tem .lt. zero .and. errq .lt. 0.1) then
3293 skpup = .true.
3294
3295 = 0.0
3296
3297
3298
3299 ELSE
3300 ERRQ = TX2
3301
3302
3303
3304 ENDIF
3305 ENDIF
3306
3307
3308
3309 ENDIF
3310
3311 ENDDO
3312
3313
3314
3315
3316
3317
3318
3319 IF (ERRQ .LT. 0.1) THEN
3320 DDFT = .TRUE.
3321 RNB = - RNB
3322
3323
3324
3325 ELSE
3326 DDFT = .FALSE.
3327 ENDIF
3328
3329
3330
3331
3332 IF (DDFT) THEN
3333 TX1 = 0.0
3334 DO L=KD,KB1
3335 TX1 = TX1 + RNF(L)
3336 ENDDO
3337
3338 = TRAIN / (TX1+RNT+RNB)
3339 IF (ABS(TX1-1.0) .LT. 0.2) THEN
3340 RNT = MAX(RNT*TX1,ZERO)
3341 RNB = RNB * TX1
3342 ELSE
3343 DDFT = .FALSE.
3344 ERRQ = 10.0
3345 ENDIF
3346 ENDIF
3347 enddo
3348
3349 = 0.0
3350 IF (.NOT. DDFT) RETURN
3351
3352
3353 DO L=KD,KB1
3354 RNF(L) = RNF(L) * TX1
3355
3356 ENDDO
3357
3358
3359
3360
3361
3362
3363
3364 DO L=KD,K
3365 WCB(L) = 0.0
3366 ENDDO
3367
3368 = .NOT. DDFT
3369
3370 = 10.0
3371 IF (.NOT. SKPDD) THEN
3372
3373
3374
3375
3376 = MAX(KB1,KD1)
3377 DO L=KK,K
3378 STLT(L) = STLT(L-1)
3379 ENDDO
3380 TEM1 = 1.0 / BB1
3381
3382 DO L=KD,K
3383 IF (L .LE. KBL) THEN
3384 TEM = STLA * TEM1
3385 STLT(L) = ETA(L) * STLT(L) * TEM / ROR(L)
3386 ELSE
3387 STLT(L) = 0.0
3388 ENDIF
3389 ENDDO
3390
3391
3392 = 0.0
3393 rsum2 = 0.0
3394
3395
3396 = 99
3397 DO L=KD,K+1
3398 ETD(L) = 0.0
3399 WVL(L) = 0.0
3400
3401 ENDDO
3402 DO L=KD,K
3403 EVP(L) = 0.0
3404 BUY(L) = 0.0
3405 QRP(L+1) = 0.0
3406 ENDDO
3407 HOD(KD) = HOL(KD)
3408 QOD(KD) = QOL(KD)
3409 TX1 = 0.0
3410
3411
3412
3413 = 0.0
3414 TX5 = TX1
3415 QA(1) = 0.0
3416
3417
3418
3419
3420
3421
3422 IF (RNT .GT. 0.0) THEN
3423 if (TX1 .gt. 0.0) THEN
3424 QRP(KD) = (RPART*RNT / (ROR(KD)*TX1*GMS(KD))) &
3425 & ** (1.0/1.1364)
3426 else
3427 tx1 = RPART*RNT / (ROR(KD)*GMS(KD)*QRP(KD)**1.1364)
3428 endif
3429 RNTP = (1.0 - RPART) * RNT
3430 BUY(KD) = - ROR(KD) * TX1 * QRP(KD)
3431 ELSE
3432 QRP(KD) = 0.0
3433 ENDIF
3434
3435
3436
3437
3438 = 1
3439 DO L=KD1,K+1
3440
3441 QA(1) = 0.0
3442 ddlgk = idn(idnm) .eq. 99
3443 if (.not. ddlgk) cycle
3444 IF (L .LE. K) THEN
3445 ST1 = 1.0 - ALFIND(L)
3446 WA(1) = ALFIND(L)*HOL(L-1) + ST1*HOL(L)
3447 WA(2) = ALFIND(L)*QOL(L-1) + ST1*QOL(L)
3448 WA(3) = ALFIND(L)*TOL(L-1) + ST1*TOL(L)
3449 QA(2) = ALFIND(L)*HST(L-1) + ST1*HST(L)
3450 QA(3) = ALFIND(L)*QST(L-1) + ST1*QST(L)
3451 ELSE
3452 WA(1) = HOL(K)
3453 WA(2) = QOL(K)
3454 WA(3) = TOL(K)
3455 QA(2) = HST(K)
3456 QA(3) = QST(K)
3457 ENDIF
3458
3459 = 2.0
3460 IF (L .EQ. KD1) FAC = 1.0
3461
3462 FACG = FAC * 0.5 * GMF5
3463
3464
3465 (KD) = ROR(L)
3466
3467
3468 = TX5
3469 WVL(L) = MAX(WVL(L-1),ONE_M1)
3470
3471 QRP(L) = MAX(QRP(L-1),QRP(L))
3472
3473
3474 (1) = GMS(L-1) * QRPF(QRP(L-1))
3475 RNT = ROR(L-1) * (WVL(L-1)+VT(1))*QRP(L-1)
3476
3477
3478
3479
3480
3481
3482
3483 = MAX(ALM,ONE_M6) * MAX(ETA(L), ONE)
3484
3485 (1) = PICON*TEM*(QRB(L-1)+QRT(L-1))
3486 TRW(2) = 1.0 / TRW(1)
3487
3488 (1) = 0.5 * (GAM(L-1) + GAM(L))
3489 VRW(2) = 1.0 / (VRW(1) + VRW(1))
3490
3491 = (QRT(L-1)+QRB(L-1))*(ONEBG*FAC*500.00*EKNOB)
3492
3493 = 1.0 / (WA(3) * (1.0 + NU*WA(2)))
3494
3495 (L) = ETD(L-1)
3496 HOD(L) = HOD(L-1)
3497 QOD(L) = QOD(L-1)
3498
3499 = 10.0
3500
3501
3502 IF (L .LE. KBL) THEN
3503 TX3 = STLT(L-1) * QRT(L-1) * (0.5*FAC)
3504 TX8 = STLT(L) * QRB(L-1) * (0.5*FAC)
3505 TX9 = TX8 + TX3
3506 ELSE
3507 TX3 = 0.0
3508 TX8 = 0.0
3509 TX9 = 0.0
3510 ENDIF
3511
3512 = WVL(L-1) + VT(1)
3513 IF (TEM .GT. 0.0) THEN
3514 TEM1 = 1.0 / (TEM*ROR(L-1))
3515 TX3 = VT(1) * TEM1 * ROR(L-1) * TX3
3516 TX6 = TX1 * TEM1
3517 ELSE
3518 TX6 = 1.0
3519 ENDIF
3520
3521
3522 IF (L .EQ. KD1) THEN
3523 IF (RNT .GT. 0.0) THEN
3524 TEM = MAX(QRP(L-1),QRP(L))
3525 WVL(L) = TX1 * TEM * QRB(L-1)*(FACG*5.0)
3526 ENDIF
3527 WVL(L) = MAX(ONE_M2, WVL(L))
3528 TRW(1) = TRW(1) * 0.5
3529 TRW(2) = TRW(2) + TRW(2)
3530 ELSE
3531 IF (DDLGK) EVP(L-1) = EVP(L-2)
3532 ENDIF
3533
3534
3535
3536
3537 IF (L .LT. IDH) THEN
3538
3539 ETD(L) = 0.0
3540 HOD(L) = WA(1)
3541 QOD(L) = WA(2)
3542 EVP(L-1) = 0.0
3543 WVL(L) = 0.0
3544 QRP(L) = 0.0
3545 BUY(L) = 0.0
3546 TX5 = TX9
3547 ERRQ = 0.0
3548 RNTP = RNTP + RNT * TX1
3549 RNT = 0.0
3550 WCB(L-1) = 0.0
3551 ENDIF
3552
3553
3554
3555
3556
3557
3558 DO ITR=1,ITRMD
3559
3560
3561 = ERRQ .GT. ERRMIN
3562 IF (cnvflg) THEN
3563
3564
3565 (1) = GMS(L) * QRPF(QRP(L))
3566 TEM = WVL(L) + VT(1)
3567
3568 IF (TEM .GT. 0.0) THEN
3569 ST1 = ROR(L) * TEM * QRP(L) + RNT
3570 IF (ST1 .NE. 0.0) ST1 = 2.0 * EVP(L-1) / ST1
3571 TEM1 = 1.0 / (TEM*ROR(L))
3572 TEM2 = VT(1) * TEM1 * ROR(L) * TX8
3573 ELSE
3574 TEM1 = 0.0
3575 TEM2 = TX8
3576 ST1 = 0.0
3577 ENDIF
3578
3579
3580
3581
3582 = tx5
3583 TEM = ROR(L)*WVL(L) - ROR(L-1)*WVL(L-1)
3584 if (tem .gt. 0.0) then
3585 TX5 = (TX1 - ST1 + TEM2 + TX3)/(1.0+tem*tem1)
3586 else
3587 TX5 = TX1 - tem*tx6 - ST1 + TEM2 + TX3
3588 endif
3589 TX5 = MAX(TX5,ZERO)
3590 tx5 = 0.5 * (tx5 + st2)
3591
3592
3593
3594
3595
3596
3597
3598
3599
3600
3601
3602
3603
3604
3605
3606
3607
3608
3609
3610
3611 = ETD(L)
3612 ETD(L) = ROR(L) * TX5 * MAX(WVL(L),ZERO)
3613
3614 if (etd(l) .gt. 0.0) etd(l) = 0.5 * (etd(l) + tem1)
3615
3616
3617 = ETD(L) - ETD(L-1)
3618
3619
3620
3621
3622
3623
3624
3625
3626
3627 = ETD(L) - TEM1
3628
3629 = max(abs(del_eta), trw(1))
3630 tem2 = del_eta / tem
3631 TEM1 = SQRT(MAX((tem+DEL_ETA)*(tem-DEL_ETA),ZERO))
3632
3633
3634 = (0.5 + ASIN(TEM2)*PIINV)*DEL_ETA + TEM1*PIINV
3635
3636 DDZ = EDZ - DEL_ETA
3637 WCB(L-1) = ETD(L) + DDZ
3638
3639 = HOD(L)
3640 IF (DEL_ETA .GT. 0.0) THEN
3641 QQQ = 1.0 / (ETD(L) + DDZ)
3642 HOD(L) = (ETD(L-1)*HOD(L-1) + DEL_ETA*HOL(L-1) &
3643 & + DDZ*WA(1)) * QQQ
3644 QOD(L) = (ETD(L-1)*QOD(L-1) + DEL_ETA*QOL(L-1) &
3645 & + DDZ*WA(2)) * QQQ
3646 ELSEif((ETD(L-1) + EDZ) .gt. 0.0) then
3647 QQQ = 1.0 / (ETD(L-1) + EDZ)
3648 HOD(L) = (ETD(L-1)*HOD(L-1) + EDZ*WA(1)) * QQQ
3649 QOD(L) = (ETD(L-1)*QOD(L-1) + EDZ*WA(2)) * QQQ
3650 ENDIF
3651 ERRH = HOD(L) - TEM1
3652 ERRQ = ABS(ERRH/HOD(L)) + ABS(ERRE/MAX(ETD(L),ONE_M5))
3653
3654
3655 = DDZ
3656 VT(2) = QQQ
3657
3658
3659 = DOF
3660 TEM4 = QOD(L)
3661 TEM1 = VRW(1)
3662
3663 = QA(3) + 0.5 * (GAF(L-1)+GAF(L)) &
3664 & * (HOD(L)-QA(2))
3665
3666
3667
3668 = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
3669 TEM2 = ROR(L) * QRP(L)
3670 CALL QRABF(TEM2,QRAF,QRBF)
3671 TEM6 = TX5 * (1.6 + 124.9 * QRAF) * QRBF * TX4
3672
3673 = TEM6 * ST2 / ((5.4E5*ST2 + 2.55E6)*(ETD(L)+DDZ))
3674
3675 = - ((1.0+TEM1)*(QHS+CE) + TEM1*QOD(L))
3676 TEM3 = (1.0 + TEM1) * QHS * (QOD(L)+CE)
3677 TEM = MAX(TEM2*TEM2 - 4.0*TEM1*TEM3,ZERO)
3678 QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
3679
3680
3681
3682
3683
3684 = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
3685 CE = TEM6 * ST2 / ((5.4E5*ST2 + 2.55E6)*(ETD(L)+DDZ))
3686
3687
3688
3689
3690 = - ((1.0+TEM1)*(QHS+CE) + TEM1*tem4)
3691 TEM3 = (1.0 + TEM1) * QHS * (tem4+CE)
3692 TEM = MAX(TEM2*TEM2 - 4.0*TEM1*TEM3,ZERO)
3693 QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
3694
3695
3696
3697 (L-1) = (QOD(L)-TEM4) * (ETD(L)+DDZ)
3698
3699 (1) = TX1*RNT + RNF(L-1) - EVP(L-1)
3700
3701
3702
3703
3704
3705 if (qa(1) .gt. 0.0) then
3706 IF (ETD(L) .GT. 0.0) THEN
3707 TEM = QA(1) / (ETD(L)+ROR(L)*TX5*VT(1))
3708 QRP(L) = MAX(TEM,ZERO)
3709 ELSEIF (TX5 .GT. 0.0) THEN
3710 QRP(L) = (MAX(ZERO,QA(1)/(ROR(L)*TX5*GMS(L)))) &
3711 & ** (1.0/1.1364)
3712 ELSE
3713 QRP(L) = 0.0
3714 ENDIF
3715 else
3716 qrp(l) = 0.5 * qrp(l)
3717 endif
3718
3719 = WA(3)+(HOD(L)-WA(1)-ALHL*(QOD(L)-WA(2))) &
3720 & * (1.0/CP)
3721
3722
3723
3724 = TEM1 * (1.0 + NU*QOD(L))
3725 ROR(L) = CMPOR * PRL(L) / TEM1
3726 TEM1 = TEM1 * DOFW
3727
3728
3729 (L) = (TEM1 - 1.0 - QRP(L)) * ROR(L) * TX5
3730
3731
3732 = WVL(L)
3733
3734 (L) = VT(2) * (ETD(L-1)*WVL(L-1) - FACG &
3735 & * (BUY(L-1)*QRT(L-1)+BUY(L)*QRB(L-1)))
3736
3737
3738
3739
3740
3741
3742 if (wvl(l) .lt. 0.0) then
3743
3744
3745
3746
3747 (L) = 1.0e-10
3748 else
3749 WVL(L) = 0.5*(WVL(L)+TEM1)
3750 endif
3751
3752
3753
3754
3755 = WVL(L) - TEM1
3756
3757 = ERRQ + ABS(ERRW/MAX(WVL(L),ONE_M5))
3758
3759
3760
3761
3762
3763
3764
3765
3766
3767
3768
3769 IF (ITR .GE. MIN(ITRMND,ITRMD/2)) THEN
3770
3771 IF (ETD(L-1) .EQ. 0.0 .AND. ERRQ .GT. 0.2) THEN
3772
3773 (L) = BUD(KD)
3774 ETD(L) = 0.0
3775 WVL(L) = 0.0
3776 ERRQ = 0.0
3777 HOD(L) = WA(1)
3778 QOD(L) = WA(2)
3779
3780 if (L .le. KBL) then
3781 TX5 = TX9
3782 else
3783 TX5 = (STLT(KB1) * QRT(KB1) &
3784 & + STLT(KBL) * QRB(KB1)) * (0.5*FAC)
3785 endif
3786
3787
3788
3789
3790 (L-1) = 0.0
3791 TEM = MAX(TX1*RNT+RNF(L-1),ZERO)
3792 QA(1) = TEM - EVP(L-1)
3793
3794
3795
3796
3797
3798
3799
3800
3801
3802
3803
3804
3805 (L) = (QA(1) / (ROR(L)*TX5*GMS(L))) &
3806 & ** (1.0/1.1364)
3807
3808 (L) = - ROR(L) * TX5 * QRP(L)
3809 WCB(L-1) = 0.0
3810 ENDIF
3811
3812 = ETD(L) - ETD(L-1)
3813 IF(DEL_ETA .LT. 0.0 .AND. ERRQ .GT. 0.1) THEN
3814 ROR(L) = BUD(KD)
3815 ETD(L) = 0.0
3816 WVL(L) = 0.0
3817
3818 (L-1) = TX5
3819
3820 = - ETD(L-1)
3821 EDZ = 0.0
3822 DDZ = -DEL_ETA
3823 WCB(L-1) = DDZ
3824
3825 (L) = HOD(L-1)
3826 QOD(L) = QOD(L-1)
3827
3828 = QOD(L)
3829 TEM1 = VRW(1)
3830
3831 = QA(3) + 0.5 * (GAF(L-1)+GAF(L)) &
3832 & * (HOD(L)-QA(2))
3833
3834
3835
3836
3837 = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
3838 TEM2 = ROR(L) * QRP(L-1)
3839 CALL QRABF(TEM2,QRAF,QRBF)
3840 TEM6 = TX5 * (1.6 + 124.9 * QRAF) * QRBF * TX4
3841
3842 = TEM6*ST2/((5.4E5*ST2 + 2.55E6)*(ETD(L)+DDZ))
3843
3844
3845 = - ((1.0+TEM1)*(QHS+CE) + TEM1*QOD(L))
3846 TEM3 = (1.0 + TEM1) * QHS * (QOD(L)+CE)
3847 TEM = MAX(TEM2*TEM2 -FOUR*TEM1*TEM3,ZERO)
3848 QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
3849
3850
3851
3852 = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
3853 CE = TEM6*ST2/((5.4E5*ST2 + 2.55E6)*(ETD(L)+DDZ))
3854
3855
3856
3857
3858 = - ((1.0+TEM1)*(QHS+CE) + TEM1*tem4)
3859 TEM3 = (1.0 + TEM1) * QHS * (tem4+CE)
3860 TEM = MAX(TEM2*TEM2 -FOUR*TEM1*TEM3,ZERO)
3861 QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
3862
3863
3864
3865 (L-1) = (QOD(L)-TEM4) * (ETD(L)+DDZ)
3866
3867
3868
3869
3870 (1) = TX1*RNT + RNF(L-1)
3871 EVP(L-1) = min(EVP(L-1), QA(1))
3872 QA(1) = QA(1) - EVP(L-1)
3873 qrp(l) = 0.0
3874
3875
3876
3877
3878
3879
3880
3881
3882
3883
3884
3885
3886
3887
3888
3889
3890
3891
3892
3893
3894
3895
3896
3897
3898 IF (L .LE. K) THEN
3899 RNS(L) = QA(1)
3900 QA(1) = 0.0
3901 ENDIF
3902 tx5 = tx9
3903 ERRQ = 0.0
3904 QRP(L) = 0.0
3905 BUY(L) = 0.0
3906
3907 ENDIF
3908 ENDIF
3909 ENDIF
3910
3911 ENDDO
3912 IF (L .LE. K) THEN
3913 IF (ETD(L-1) .EQ. 0.0 &
3914 & .AND. ERRQ .GT. 0.1 .and. l .le. kbl) THEN
3915
3916
3917 (L) = BUD(KD)
3918 HOD(L) = WA(1)
3919 QOD(L) = WA(2)
3920 TX5 = TX9
3921
3922 (L-1) = 0.0
3923
3924 (1) = TX1*RNT + RNF(L-1)
3925 EVP(L-1) = min(EVP(L-1), QA(1))
3926 QA(1) = QA(1) - EVP(L-1)
3927
3928
3929
3930
3931
3932
3933
3934
3935
3936 (L) = (QA(1) / (ROR(L)*TX5*GMS(L))) &
3937 & ** (1.0/1.1364)
3938
3939 (L) = 0.0
3940 WVL(L) = 0.0
3941 ST1 = 1.0 - ALFIND(L)
3942
3943 ERRQ = 0.0
3944 BUY(L) = - ROR(L) * TX5 * QRP(L)
3945 WCB(L-1) = 0.0
3946 ENDIF
3947 ENDIF
3948
3949 = MIN(IDN(idnm), K+1)
3950 IF (ERRQ .LT. 1.0 .AND. L .LE. LL) THEN
3951 IF (ETD(L-1) .GT. 0.0 .AND. ETD(L) .EQ. 0.0) THEN
3952 IDN(idnm) = L
3953 wvl(l) = 0.0
3954 if (L .lt. KBL .or. tx5 .gt. 0.0) idnm = idnm + 1
3955 errq = 0.0
3956 ENDIF
3957 if (etd(l) .eq. 0.0 .and. l .gt. kbl) then
3958 idn(idnm) = l
3959 if (tx5 .gt. 0.0) idnm = idnm + 1
3960 endif
3961 ENDIF
3962
3963
3964
3965
3966
3967
3968
3969
3970
3971
3972
3973
3974 IF (ERRQ .GT. 0.1 .AND. IDN(idnm) .EQ. 99) &
3975 & DDFT = .FALSE.
3976
3977
3978 = 0.0
3979 IF (.NOT. DDFT) RETURN
3980
3981
3982
3983
3984
3985
3986
3987
3988
3989
3990
3991 ENDDO
3992
3993 = 0.0
3994
3995 DOF = QA(1)
3996
3997
3998
3999
4000 ENDIF
4001
4002
4003 (KD) = RNTP
4004 TX1 = EVP(KD)
4005 TX2 = RNTP + RNB + DOF
4006
4007
4008 = IDH
4009 IF (II .GE. KD1+1) THEN
4010 RNN(KD) = RNN(KD) + RNF(KD)
4011 TX2 = TX2 + RNF(KD)
4012 RNN(II-1) = 0.0
4013 TX1 = EVP(II-1)
4014 ENDIF
4015
4016 DO L=KD,K
4017 II = IDH
4018
4019 IF (L .GT. KD1 .AND. L .LT. II) THEN
4020 RNN(L-1) = RNF(L-1)
4021 TX2 = TX2 + RNN(L-1)
4022 ELSEIF (L .GE. II .AND. L .LT. IDN(idnm)) THEN
4023 rnn(l) = rns(l)
4024 tx2 = tx2 + rnn(l)
4025 TX1 = TX1 + EVP(L)
4026 ELSEIF (L .GE. IDN(idnm)) THEN
4027 ETD(L+1) = 0.0
4028 HOD(L+1) = 0.0
4029 QOD(L+1) = 0.0
4030 EVP(L) = 0.0
4031 RNN(L) = RNF(L) + RNS(L)
4032 TX2 = TX2 + RNN(L)
4033 ENDIF
4034
4035 ENDDO
4036
4037
4038
4039 = KBL
4040
4041 RNN(L) = RNN(L) + RNB
4042 CLDFRD(L) = TX5
4043
4044
4045
4046
4047
4048
4049
4050
4051 IF (TX1 .GT. 0.0) THEN
4052 TX1 = (TRAIN - TX2) / TX1
4053 ELSE
4054 TX1 = 0.0
4055 ENDIF
4056
4057 DO L=KD,K
4058 EVP(L) = EVP(L) * TX1
4059 ENDDO
4060
4061
4062
4063
4064 RETURN
4065 END
4066
4067 SUBROUTINE QSATCN(TT,P,Q,DQDT)
4068
4069
4070 USE MACHINE , ONLY : kind_phys
4071 USE FUNCPHYS , ONLY : fpvs
4072 USE PHYSCONS, RV => con_RV, CVAP => con_CVAP, CLIQ => con_CLIQ &
4073 &, CSOL => con_CSOL, TTP => con_TTP, HVAP => con_HVAP &
4074 &, HFUS => con_HFUS, EPS => con_eps, EPSM1 => con_epsm1
4075 implicit none
4076
4077 real(kind=kind_phys) TT, P, Q, DQDT
4078
4079 real(kind=kind_phys) rvi, facw, faci, hsub, tmix, DEN
4080 real(kind=kind_phys) ZERO,ONE,ONE_M10
4081 PARAMETER (RVI=1.0/RV)
4082 PARAMETER (FACW=CVAP-CLIQ, FACI=CVAP-CSOL)
4083 PARAMETER (HSUB=HVAP+HFUS, tmix=TTP-20.0, DEN=1.0/(TTP-TMIX))
4084 PARAMETER (ZERO=0.,ONE=1.,ONE_M10=1.E-10)
4085
4086
4087 real(kind=kind_phys) es, d, hlorv, W
4088
4089
4090 = 0.01 * fpvs(tt)
4091 = 1.0 / max(p+epsm1*es,ONE_M10)
4092
4093 = MIN(eps*es*D, ONE)
4094
4095 = max(ZERO, min(ONE, (TT - TMIX)*DEN))
4096 hlorv = ( W * (HVAP + FACW * (tt-ttp)) &
4097 & + (1.0-W) * (HSUB + FACI * (tt-ttp)) ) * RVI
4098 dqdt = p * q * hlorv * D / (tt*tt)
4099
4100 return
4101 end
4102
4103 SUBROUTINE ANGRAD( PRES, ALM, AL2, TLA, PRB, WFN, UFN)
4104 USE MACHINE , ONLY : kind_phys
4105 use module_ras , only : refp, refr, tlac, plac, tlbpl, drdp, almax
4106 implicit none
4107
4108 real(kind=kind_phys) PRES &
4109 &, ALM, AL2, TLA, TEM, TEM1 &
4110 &, PRB, ACR, WFN, UFN
4111
4112 integer i
4113
4114 IF (TLA .LT. 0.0) THEN
4115 IF (PRES .LE. PLAC(1)) THEN
4116 TLA = TLAC(1)
4117 ELSEIF (PRES .LE. PLAC(2)) THEN
4118 TLA = TLAC(2) + (PRES-PLAC(2))*tlbpl(1)
4119 ELSEIF (PRES .LE. PLAC(3)) THEN
4120 TLA = TLAC(3) + (PRES-PLAC(3))*tlbpl(2)
4121 ELSEIF (PRES .LE. PLAC(4)) THEN
4122 TLA = TLAC(4) + (PRES-PLAC(4))*tlbpl(3)
4123 ELSEIF (PRES .LE. PLAC(5)) THEN
4124 TLA = TLAC(5) + (PRES-PLAC(5))*tlbpl(4)
4125 ELSEIF (PRES .LE. PLAC(6)) THEN
4126 TLA = TLAC(6) + (PRES-PLAC(6))*tlbpl(5)
4127 ELSEIF (PRES .LE. PLAC(7)) THEN
4128 TLA = TLAC(7) + (PRES-PLAC(7))*tlbpl(6)
4129 ELSEIF (PRES .LE. PLAC(8)) THEN
4130 TLA = TLAC(8) + (PRES-PLAC(8))*tlbpl(7)
4131 ELSE
4132 TLA = TLAC(8)
4133 ENDIF
4134 ENDIF
4135 IF (PRES .GE. REFP(1)) THEN
4136 TEM = REFR(1)
4137 ELSEIF (PRES .GE. REFP(2)) THEN
4138 TEM = REFR(1) + (PRES-REFP(1)) * drdp(1)
4139 ELSEIF (PRES .GE. REFP(3)) THEN
4140 TEM = REFR(2) + (PRES-REFP(2)) * drdp(2)
4141 ELSEIF (PRES .GE. REFP(4)) THEN
4142 TEM = REFR(3) + (PRES-REFP(3)) * drdp(3)
4143 ELSEIF (PRES .GE. REFP(5)) THEN
4144 TEM = REFR(4) + (PRES-REFP(4)) * drdp(4)
4145 ELSEIF (PRES .GE. REFP(6)) THEN
4146 TEM = REFR(5) + (PRES-REFP(5)) * drdp(5)
4147 ELSE
4148 TEM = REFR(6)
4149 ENDIF
4150
4151 = 2.0E-4 / tem
4152 al2 = min(4.0*tem, max(alm, tem))
4153
4154 RETURN
4155 END
4156 SUBROUTINE SETQRP
4157 USE MACHINE , ONLY : kind_phys
4158 use module_ras , only : NQRP,C1XQRP,C2XQRP,TBQRP,TBQRA,TBQRB
4159 implicit none
4160
4161 real(kind=kind_phys) tem2,tem1,x,xinc,xmax,xmin
4162 integer jx
4163
4164
4165 = 0.0
4166 XMAX = 5.0
4167 XINC = (XMAX-XMIN)/(NQRP-1)
4168 C2XQRP = 1.0/XINC
4169 C1XQRP = 1.0 - XMIN*C2XQRP
4170 TEM1 = 0.001 ** 0.2046
4171 TEM2 = 0.001 ** 0.525
4172 DO JX=1,NQRP
4173 X = XMIN + (JX-1)*XINC
4174 TBQRP(JX) = X ** 0.1364
4175 TBQRA(JX) = TEM1 * X ** 0.2046
4176 TBQRB(JX) = TEM2 * X ** 0.525
4177 ENDDO
4178
4179 RETURN
4180 END
4181 FUNCTION QRPF(QRP)
4182
4183 USE MACHINE , ONLY : kind_phys
4184 use module_ras , only : NQRP,C1XQRP,C2XQRP,TBQRP,TBQRA,TBQRB
4185 implicit none
4186
4187 real(kind=kind_phys) QRP, QRPF, XJ, REAL_NQRP, ONE
4188 PARAMETER (ONE=1.0)
4189 INTEGER JX
4190
4191 = REAL(NQRP)
4192 XJ = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),REAL_NQRP)
4193
4194 = MIN(XJ,NQRP-ONE)
4195 QRPF = TBQRP(JX) + (XJ-JX) * (TBQRP(JX+1)-TBQRP(JX))
4196
4197 RETURN
4198 END
4199 SUBROUTINE QRABF(QRP,QRAF,QRBF)
4200 USE MACHINE , ONLY : kind_phys
4201 use module_ras , only : NQRP,C1XQRP,C2XQRP,TBQRP,TBQRA,TBQRB
4202 implicit none
4203
4204 real(kind=kind_phys) QRP, QRAF, QRBF, XJ, REAL_NQRP, ONE
4205 PARAMETER (ONE=1.0)
4206 INTEGER JX
4207
4208 =REAL(NQRP)
4209 XJ = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),REAL_NQRP)
4210 JX = MIN(XJ,NQRP-ONE)
4211 XJ = XJ - JX
4212 QRAF = TBQRA(JX) + XJ * (TBQRA(JX+1)-TBQRA(JX))
4213 QRBF = TBQRB(JX) + XJ * (TBQRB(JX+1)-TBQRB(JX))
4214
4215 RETURN
4216 END
4217 SUBROUTINE SETVTP
4218 USE MACHINE , ONLY : kind_phys
4219 use module_ras , only : NVTP,C1XVTP,C2XVTP,TBVTP
4220 implicit none
4221
4222 real(kind=kind_phys) vtpexp,xinc,x,xmax,xmin
4223 integer jx
4224 PARAMETER(VTPEXP=-0.3636)
4225
4226 = 0.05
4227 XMAX = 1.5
4228 XINC = (XMAX-XMIN)/(NVTP-1)
4229 C2XVTP = 1.0/XINC
4230 C1XVTP = 1.0 - XMIN*C2XVTP
4231 DO JX=1,NVTP
4232 X = XMIN + (JX-1)*XINC
4233 TBVTP(JX) = X ** VTPEXP
4234 ENDDO
4235
4236 RETURN
4237 END
4238 FUNCTION VTPF(ROR)
4239
4240 USE MACHINE , ONLY : kind_phys
4241 use module_ras , only : NVTP,C1XVTP,C2XVTP,TBVTP
4242 implicit none
4243 real(kind=kind_phys) ROR, VTPF, XJ, REAL_NVTP, ONE
4244 PARAMETER (ONE=1.0)
4245 INTEGER JX
4246
4247 = REAL(NVTP)
4248 XJ = MIN(MAX(C1XVTP+C2XVTP*ROR,ONE),REAL_NVTP)
4249 JX = MIN(XJ,NVTP-ONE)
4250 VTPF = TBVTP(JX) + (XJ-JX) * (TBVTP(JX+1)-TBVTP(JX))
4251
4252 RETURN
4253 END
4254 FUNCTION CLF(PRATE)
4255
4256 USE MACHINE , ONLY : kind_phys
4257 implicit none
4258 real(kind=kind_phys) PRATE, CLF
4259
4260 real (kind=kind_phys), parameter :: ccf1=0.30, ccf2=0.09 &
4261 &, ccf3=0.04, ccf4=0.01 &
4262 &, pr1=1.0, pr2=5.0 &
4263 &, pr3=20.0
4264
4265 if (prate .lt. pr1) then
4266 clf = ccf1
4267 elseif (prate .lt. pr2) then
4268 clf = ccf2
4269 elseif (prate .lt. pr3) then
4270 clf = ccf3
4271 else
4272 clf = ccf4
4273 endif
4274
4275 RETURN
4276 END
4277