1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
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
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
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
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
|
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
!**************************************************************************
!**** SUBROUTINE CONVECT *****
!**** VERSION 4.3c *****
!**** 20 May, 2002 *****
!**** Kerry Emanuel *****
!**************************************************************************
!
SUBROUTINE CONVECT &
(ND, NL, DELT, IFLAG, &
PRECIP, WD, TPRIME, QPRIME, CBMF )
!
!-cv *************************************************************************
!-cv C. Forster, November 2003 - May 2004:
!-cv
!-cv The subroutine has been downloaded from Kerry Emanuel's homepage,
!-cv where further infos on the convection scheme can be found
!-cv http://www-paoc.mit.edu/~emanuel/home.html
!-cv
!-cv The following changes have been made to integrate this subroutine
!-cv into FLEXPART
!-cv
!-cv Putting most of the variables in a new common block
!-cv renaming eps to eps0 because there is some eps already in includepar
!-cv
!-cv removing the arrays U,V,TRA and related arrays
!-cv
!-cv renaming the original arrays T,Q,QS,P,PH to
!-cv TCONV,QCONV,QSCONV,PCONV_HPA,PHCONV_HPA
!-cv
!-cv Initialization of variables has been put into parameter statements
!-cv instead of assignment of values at each call, in order to save
!-cv computation time.
!***************************************************************************
!
!-----------------------------------------------------------------------------
! *** On input: ***
!
!T: Array of absolute temperature (K) of dimension ND, with first
! index corresponding to lowest model level. Note that this array
! will be altered by the subroutine if dry convective adjustment
! occurs and if IPBL is not equal to 0.
!
!Q: Array of specific humidity (gm/gm) of dimension ND, with first
! index corresponding to lowest model level. Must be defined
! at same grid levels as T. Note that this array will be altered
! if dry convective adjustment occurs and if IPBL is not equal to 0.
!
!QS: Array of saturation specific humidity of dimension ND, with first
! index corresponding to lowest model level. Must be defined
! at same grid levels as T. Note that this array will be altered
! if dry convective adjustment occurs and if IPBL is not equal to 0.
!
!U: Array of zonal wind velocity (m/s) of dimension ND, witth first
! index corresponding with the lowest model level. Defined at
! same levels as T. Note that this array will be altered if
! dry convective adjustment occurs and if IPBL is not equal to 0.
!
!V: Same as U but for meridional velocity.
!
!TRA: Array of passive tracer mixing ratio, of dimensions (ND,NTRA),
! where NTRA is the number of different tracers. If no
! convective tracer transport is needed, define a dummy
! input array of dimension (ND,1). Tracers are defined at
! same vertical levels as T. Note that this array will be altered
! if dry convective adjustment occurs and if IPBL is not equal to 0.
!
!P: Array of pressure (mb) of dimension ND, with first
! index corresponding to lowest model level. Must be defined
! at same grid levels as T.
!
!PH: Array of pressure (mb) of dimension ND+1, with first index
! corresponding to lowest level. These pressures are defined at
! levels intermediate between those of P, T, Q and QS. The first
! value of PH should be greater than (i.e. at a lower level than)
! the first value of the array P.
!
!ND: The dimension of the arrays T,Q,QS,P,PH,FT and FQ
!
!NL: The maximum number of levels to which convection can
! penetrate, plus 1.
! NL MUST be less than or equal to ND-1.
!
!NTRA:The number of different tracers. If no tracer transport
! is needed, set this equal to 1. (On most compilers, setting
! NTRA to 0 will bypass tracer calculation, saving some CPU.)
!
!DELT: The model time step (sec) between calls to CONVECT
!
!----------------------------------------------------------------------------
! *** On Output: ***
!
!IFLAG: An output integer whose value denotes the following:
!
! VALUE INTERPRETATION
! ----- --------------
! 0 No moist convection; atmosphere is not
! unstable, or surface temperature is less
! than 250 K or surface specific humidity
! is non-positive.
!
! 1 Moist convection occurs.
!
! 2 No moist convection: lifted condensation
! level is above the 200 mb level.
!
! 3 No moist convection: cloud base is higher
! then the level NL-1.
!
! 4 Moist convection occurs, but a CFL condition
! on the subsidence warming is violated. This
! does not cause the scheme to terminate.
!
!FT: Array of temperature tendency (K/s) of dimension ND, defined at same
! grid levels as T, Q, QS and P.
!
!FQ: Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,
! defined at same grid levels as T, Q, QS and P.
!
!FU: Array of forcing of zonal velocity (m/s^2) of dimension ND,
! defined at same grid levels as T.
!
!FV: Same as FU, but for forcing of meridional velocity.
!
!FTRA: Array of forcing of tracer content, in tracer mixing ratio per
! second, defined at same levels as T. Dimensioned (ND,NTRA).
!
!PRECIP: Scalar convective precipitation rate (mm/day).
!
!WD: A convective downdraft velocity scale. For use in surface
! flux parameterizations. See convect.ps file for details.
!
!TPRIME: A convective downdraft temperature perturbation scale (K).
! For use in surface flux parameterizations. See convect.ps
! file for details.
!
!QPRIME: A convective downdraft specific humidity
! perturbation scale (gm/gm).
! For use in surface flux parameterizations. See convect.ps
! file for details.
!
!CBMF: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST
! BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT
! ITS NEXT CALL. That is, the value of CBMF must be "remembered"
! by the calling program between calls to CONVECT.
!
!-----------------------------------------------------------------------------
!
! *** THE PARAMETER NA SHOULD IN GENERAL BE GREATER THAN ***
! *** OR EQUAL TO ND + 1 ***
!
!
use par_mod
use conv_mod
implicit none
!
!-cv====>Begin Module CONVECT File convect.f Undeclared variables
!
!Argument variables
!
integer :: iflag, nd, nl
!
real :: cbmf, delt, precip, qprime, tprime, wd
!
!Local variables
!
integer :: i, icb, ihmin, inb, inb1, j, jtt, k
integer :: nk
!
real :: ad, afac, ahmax, ahmin, alt, altem
real :: am, amp1, anum, asij, awat, b6, bf2, bsum, by
real :: byp, c6, cape, capem, cbmfold, chi, coeff
real :: cpinv, cwat, damps, dbo, dbosum
real :: defrac, dei, delm, delp, delt0, delti, denom, dhdp
real :: dpinv, dtma, dtmin, dtpbl, elacrit, ents
real :: epmax, fac, fqold, frac, ftold
real :: plcl, qp1, qsm, qstm, qti, rat
real :: rdcp, revap, rh, scrit, sigt, sjmax
real :: sjmin, smid, smin, stemp, tca
real :: tvaplcl, tvpplcl, tvx, tvy, wdtrain
!integer jc,jn
!real alvnew,a2,ahm,alv,rm,sum,qnew,dphinv,tc,thbar,tnew,x
real :: FUP(NA),FDOWN(NA)
!
!-cv====>End Module CONVECT File convect.f
INTEGER :: NENT(NA)
REAL :: M(NA),MP(NA),MENT(NA,NA),QENT(NA,NA),ELIJ(NA,NA)
REAL :: SIJ(NA,NA),TVP(NA),TV(NA),WATER(NA)
REAL :: QP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA)
REAL :: SIGP(NA),TP(NA),CPN(NA)
REAL :: LV(NA),LVCP(NA),H(NA),HP(NA),GZ(NA),HM(NA)
!REAL TOLD(NA)
!
! -----------------------------------------------------------------------
!
! *** Specify Switches ***
!
! *** IPBL: Set to zero to bypass dry adiabatic adjustment ***
! *** Any other value results in dry adiabatic adjustment ***
! *** (Zero value recommended for use in models with ***
! *** boundary layer schemes) ***
!
! *** MINORIG: Lowest level from which convection may originate ***
! *** (Should be first model level at which T is defined ***
! *** for models using bulk PBL schemes; otherwise, it should ***
! *** be the first model level at which T is defined above ***
! *** the surface layer) ***
!
INTEGER,PARAMETER :: IPBL=0
INTEGER,PARAMETER :: MINORIG=1
!
!------------------------------------------------------------------------------
!
! *** SPECIFY PARAMETERS ***
!
! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
! *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- ***
! *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO ***
! *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY ***
! *** BETWEEN 0 C AND TLCRIT) ***
! *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT ***
! *** FORMULATION ***
! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT ***
! *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
! *** OF CLOUD ***
! *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN ***
! *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW ***
! *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION ***
! *** OF RAIN ***
! *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION ***
! *** OF SNOW ***
! *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM ***
! *** TRANSPORT ***
! *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION ***
! *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC ***
! *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF ***
! *** APPROACH TO QUASI-EQUILIBRIUM ***
! *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) ***
! *** (DAMP MUST BE LESS THAN 1) ***
!
REAL,PARAMETER :: ELCRIT=.0011
REAL,PARAMETER :: TLCRIT=-55.0
REAL,PARAMETER :: ENTP=1.5
REAL,PARAMETER :: SIGD=0.05
REAL,PARAMETER :: SIGS=0.12
REAL,PARAMETER :: OMTRAIN=50.0
REAL,PARAMETER :: OMTSNOW=5.5
REAL,PARAMETER :: COEFFR=1.0
REAL,PARAMETER :: COEFFS=0.8
REAL,PARAMETER :: CU=0.7
REAL,PARAMETER :: BETA=10.0
REAL,PARAMETER :: DTMAX=0.9
REAL,PARAMETER :: ALPHA=0.025 !original 0.2
REAL,PARAMETER :: DAMP=0.1
!
! *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS, ***
! *** GRAVITY, AND LIQUID WATER DENSITY. ***
! *** THESE SHOULD BE CONSISTENT WITH ***
! *** THOSE USED IN CALLING PROGRAM ***
! *** NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT ***
!
REAL,PARAMETER :: CPD=1005.7
REAL,PARAMETER :: CPV=1870.0
REAL,PARAMETER :: CL=2500.0
REAL,PARAMETER :: RV=461.5
REAL,PARAMETER :: RD=287.04
REAL,PARAMETER :: LV0=2.501E6
REAL,PARAMETER :: G=9.81
REAL,PARAMETER :: ROWL=1000.0
!
REAL,PARAMETER :: CPVMCL=CL-CPV
REAL,PARAMETER :: EPS0=RD/RV
REAL,PARAMETER :: EPSI=1./EPS0
REAL,PARAMETER :: GINV=1.0/G
REAL,PARAMETER :: EPSILON=1.e-20
! EPSILON IS A SMALL NUMBER USED TO EXCLUDE MASS FLUXES OF ZERO
!
DELTI=1.0/DELT
!
! *** INITIALIZE OUTPUT ARRAYS AND PARAMETERS ***
!
DO I=1,NL+1
FT(I)=0.0
FQ(I)=0.0
FDOWN(I)=0.0
SUB(I)=0.0
FUP(I)=0.0
M(I)=0.0
MP(I)=0.0
DO J=1,NL+1
FMASS(I,J)=0.0
MENT(I,J)=0.0
END DO
END DO
DO I=1,NL+1
RDCP=(RD*(1.-QCONV(I))+QCONV(I)*RV)/ &
(CPD*(1.-QCONV(I))+QCONV(I)*CPV)
TH(I)=TCONV(I)*(1000.0/PCONV_HPA(I))**RDCP
END DO
PRECIP=0.0
WD=0.0
TPRIME=0.0
QPRIME=0.0
IFLAG=0
!
! IF(IPBL.NE.0)THEN
!
!*** PERFORM DRY ADIABATIC ADJUSTMENT ***
!
! JC=0
! DO 30 I=NL-1,1,-1
! JN=0
! SUM=TH(I)*(1.+QCONV(I)*EPSI-QCONV(I))
! DO 10 J=I+1,NL
! SUM=SUM+TH(J)*(1.+QCONV(J)*EPSI-QCONV(J))
! THBAR=SUM/REAL(J+1-I)
! IF((TH(J)*(1.+QCONV(J)*EPSI-QCONV(J))).LT.THBAR)JN=J
! 10 CONTINUE
! IF(I.EQ.1)JN=MAX(JN,2)
! IF(JN.EQ.0)GOTO 30
! 12 CONTINUE
! AHM=0.0
! RM=0.0
! DO 15 J=I,JN
! AHM=AHM+(CPD*(1.-QCONV(J))+QCONV(J)*CPV)*TCONV(J)*
! + (PHCONV_HPA(J)-PHCONV_HPA(J+1))
! RM=RM+QCONV(J)*(PHCONV_HPA(J)-PHCONV_HPA(J+1))
! 15 CONTINUE
! DPHINV=1./(PHCONV_HPA(I)-PHCONV_HPA(JN+1))
! RM=RM*DPHINV
! A2=0.0
! DO 20 J=I,JN
! QCONV(J)=RM
! RDCP=(RD*(1.-QCONV(J))+QCONV(J)*RV)/
! 1 (CPD*(1.-QCONV(J))+QCONV(J)*CPV)
! X=(0.001*PCONV_HPA(J))**RDCP
! TOLD(J)=TCONV(J)
! TCONV(J)=X
! A2=A2+(CPD*(1.-QCONV(J))+QCONV(J)*CPV)*X*
! 1 (PHCONV_HPA(J)-PHCONV_HPA(J+1))
! 20 CONTINUE
! DO 25 J=I,JN
! TH(J)=AHM/A2
! TCONV(J)=TCONV(J)*TH(J)
! TC=TOLD(J)-273.15
! ALV=LV0-CPVMCL*TC
! QSCONV(J)=QSCONV(J)+QSCONV(J)*(1.+QSCONV(J)*(EPSI-1.))*ALV*
! 1 (TCONV(J)- TOLD(J))/(RV*TOLD(J)*TOLD(J))
! if (qslev(j) .lt. 0.) then
! write(*,*) 'qslev.lt.0 ',j,qslev
! endif
! 25 CONTINUE
! IF((TH(JN+1)*(1.+QCONV(JN+1)*EPSI-QCONV(JN+1))).LT.
! 1 (TH(JN)*(1.+QCONV(JN)*EPSI-QCONV(JN))))THEN
! JN=JN+1
! GOTO 12
! END IF
! IF(I.EQ.1)JC=JN
! 30 CONTINUE
!
! *** Remove any supersaturation that results from adjustment ***
!
!IF(JC.GT.1)THEN
! DO 38 J=1,JC
! IF(QSCONV(J).LT.QCONV(J))THEN
! ALV=LV0-CPVMCL*(TCONV(J)-273.15)
! TNEW=TCONV(J)+ALV*(QCONV(J)-QSCONV(J))/(CPD*(1.-QCONV(J))+
! 1 CL*QCONV(J)+QSCONV(J)*(CPV-CL+ALV*ALV/(RV*TCONV(J)*TCONV(J))))
! ALVNEW=LV0-CPVMCL*(TNEW-273.15)
! QNEW=(ALV*QCONV(J)-(TNEW-TCONV(J))*(CPD*(1.-QCONV(J))
! 1 +CL*QCONV(J)))/ALVNEW
! PRECIP=PRECIP+24.*3600.*1.0E5*(PHCONV_HPA(J)-PHCONV_HPA(J+1))*
! 1 (QCONV(J)-QNEW)/(G*DELT*ROWL)
! TCONV(J)=TNEW
! QCONV(J)=QNEW
! QSCONV(J)=QNEW
! END IF
! 38 CONTINUE
!END IF
!
!END IF
!
! *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY
!
GZ(1)=0.0
CPN(1)=CPD*(1.-QCONV(1))+QCONV(1)*CPV
H(1)=TCONV(1)*CPN(1)
LV(1)=LV0-CPVMCL*(TCONV(1)-273.15)
HM(1)=LV(1)*QCONV(1)
TV(1)=TCONV(1)*(1.+QCONV(1)*EPSI-QCONV(1))
AHMIN=1.0E12
IHMIN=NL
DO I=2,NL+1
TVX=TCONV(I)*(1.+QCONV(I)*EPSI-QCONV(I))
TVY=TCONV(I-1)*(1.+QCONV(I-1)*EPSI-QCONV(I-1))
GZ(I)=GZ(I-1)+0.5*RD*(TVX+TVY)*(PCONV_HPA(I-1)-PCONV_HPA(I))/ &
PHCONV_HPA(I)
CPN(I)=CPD*(1.-QCONV(I))+CPV*QCONV(I)
H(I)=TCONV(I)*CPN(I)+GZ(I)
LV(I)=LV0-CPVMCL*(TCONV(I)-273.15)
HM(I)=(CPD*(1.-QCONV(I))+CL*QCONV(I))*(TCONV(I)-TCONV(1))+ &
LV(I)*QCONV(I)+GZ(I)
TV(I)=TCONV(I)*(1.+QCONV(I)*EPSI-QCONV(I))
!
! *** Find level of minimum moist static energy ***
!
IF(I.GE.MINORIG.AND.HM(I).LT.AHMIN.AND.HM(I).LT.HM(I-1))THEN
AHMIN=HM(I)
IHMIN=I
END IF
END DO
IHMIN=MIN(IHMIN, NL-1)
!
! *** Find that model level below the level of minimum moist ***
! *** static energy that has the maximum value of moist static energy ***
!
AHMAX=0.0
! *** bug fixed: need to assign an initial value to NK
! HSO, 05.08.2009
NK=MINORIG
DO I=MINORIG,IHMIN
IF(HM(I).GT.AHMAX)THEN
NK=I
AHMAX=HM(I)
END IF
END DO
!
! *** CHECK WHETHER PARCEL LEVEL TEMPERATURE AND SPECIFIC HUMIDITY ***
! *** ARE REASONABLE ***
! *** Skip convection if HM increases monotonically upward ***
!
IF(TCONV(NK).LT.250.0.OR.QCONV(NK).LE.0.0.OR.IHMIN.EQ.(NL-1)) &
THEN
IFLAG=0
CBMF=0.0
RETURN
END IF
!
! *** CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT PARCEL ORIGIN LEVEL ***
! *** (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980) ***
!
RH=QCONV(NK)/QSCONV(NK)
CHI=TCONV(NK)/(1669.0-122.0*RH-TCONV(NK))
PLCL=PCONV_HPA(NK)*(RH**CHI)
IF(PLCL.LT.200.0.OR.PLCL.GE.2000.0)THEN
IFLAG=2
CBMF=0.0
RETURN
END IF
!
! *** CALCULATE FIRST LEVEL ABOVE LCL (=ICB) ***
!
ICB=NL-1
DO I=NK+1,NL
IF(PCONV_HPA(I).LT.PLCL)THEN
ICB=MIN(ICB,I)
END IF
END DO
IF(ICB.GE.(NL-1))THEN
IFLAG=3
CBMF=0.0
RETURN
END IF
!
! *** FIND TEMPERATURE UP THROUGH ICB AND TEST FOR INSTABILITY ***
!
! *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL ***
! *** TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC ***
! *** LIQUID WATER CONTENT ***
!
CALL TLIFT(GZ,ICB,NK,TVP,TP,CLW,ND,NL,1)
DO I=NK,ICB
TVP(I)=TVP(I)-TP(I)*QCONV(NK)
END DO
!
! *** If there was no convection at last time step and parcel ***
! *** is stable at ICB then skip rest of calculation ***
!
IF(CBMF.EQ.0.0.AND.TVP(ICB).LE.(TV(ICB)-DTMAX))THEN
IFLAG=0
RETURN
END IF
!
! *** IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY ***
!
IF(IFLAG.NE.4)IFLAG=1
!
! *** FIND THE REST OF THE LIFTED PARCEL TEMPERATURES ***
!
CALL TLIFT(GZ,ICB,NK,TVP,TP,CLW,ND,NL,2)
!
! *** SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF ***
! *** PRECIPITATION FALLING OUTSIDE OF CLOUD ***
! *** THESE MAY BE FUNCTIONS OF TP(I), PCONV_HPA(I) AND CLW(I) ***
!
DO I=1,NK
EP(I)=0.0
SIGP(I)=SIGS
END DO
DO I=NK+1,NL
TCA=TP(I)-273.15
IF(TCA.GE.0.0)THEN
ELACRIT=ELCRIT
ELSE
ELACRIT=ELCRIT*(1.0-TCA/TLCRIT)
END IF
ELACRIT=MAX(ELACRIT,0.0)
EPMAX=0.999
EP(I)=EPMAX*(1.0-ELACRIT/MAX(CLW(I),1.0E-8))
EP(I)=MAX(EP(I),0.0)
EP(I)=MIN(EP(I),EPMAX)
SIGP(I)=SIGS
END DO
!
! *** CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL ***
! *** VIRTUAL TEMPERATURE ***
!
DO I=ICB+1,NL
TVP(I)=TVP(I)-TP(I)*QCONV(NK)
END DO
TVP(NL+1)=TVP(NL)-(GZ(NL+1)-GZ(NL))/CPD
!
! *** NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS ***
!
DO I=1,NL+1
HP(I)=H(I)
NENT(I)=0
WATER(I)=0.0
EVAP(I)=0.0
WT(I)=OMTSNOW
LVCP(I)=LV(I)/CPN(I)
DO J=1,NL+1
QENT(I,J)=QCONV(J)
ELIJ(I,J)=0.0
SIJ(I,J)=0.0
END DO
END DO
QP(1)=QCONV(1)
DO I=2,NL+1
QP(I)=QCONV(I-1)
END DO
!
! *** FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S ***
! *** HIGHEST LEVEL OF NEUTRAL BUOYANCY ***
! *** AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB) ***
!
CAPE=0.0
CAPEM=0.0
INB=ICB+1
INB1=INB
BYP=0.0
DO I=ICB+1,NL-1
BY=(TVP(I)-TV(I))*(PHCONV_HPA(I)-PHCONV_HPA(I+1))/PCONV_HPA(I)
CAPE=CAPE+BY
IF(BY.GE.0.0)INB1=I+1
IF(CAPE.GT.0.0)THEN
INB=I+1
BYP=(TVP(I+1)-TV(I+1))*(PHCONV_HPA(I+1)-PHCONV_HPA(I+2))/ &
PCONV_HPA(I+1)
CAPEM=CAPE
END IF
END DO
INB=MAX(INB,INB1)
CAPE=CAPEM+BYP
DEFRAC=CAPEM-CAPE
DEFRAC=MAX(DEFRAC,0.001)
FRAC=-CAPE/DEFRAC
FRAC=MIN(FRAC,1.0)
FRAC=MAX(FRAC,0.0)
!
! *** CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL ***
!
DO I=ICB,INB
HP(I)=H(NK)+(LV(I)+(CPD-CPV)*TCONV(I))*EP(I)*CLW(I)
END DO
!
! *** CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I), ***
! *** AT EACH MODEL LEVEL ***
!
DBOSUM=0.0
!
! *** INTERPOLATE DIFFERENCE BETWEEN LIFTED PARCEL AND ***
! *** ENVIRONMENTAL TEMPERATURES TO LIFTED CONDENSATION LEVEL ***
!
TVPPLCL=TVP(ICB-1)-RD*TVP(ICB-1)*(PCONV_HPA(ICB-1)-PLCL)/ &
(CPN(ICB-1)*PCONV_HPA(ICB-1))
TVAPLCL=TV(ICB)+(TVP(ICB)-TVP(ICB+1))*(PLCL-PCONV_HPA(ICB))/ &
(PCONV_HPA(ICB)-PCONV_HPA(ICB+1))
DTPBL=0.0
DO I=NK,ICB-1
DTPBL=DTPBL+(TVP(I)-TV(I))*(PHCONV_HPA(I)-PHCONV_HPA(I+1))
END DO
DTPBL=DTPBL/(PHCONV_HPA(NK)-PHCONV_HPA(ICB))
DTMIN=TVPPLCL-TVAPLCL+DTMAX+DTPBL
DTMA=DTMIN
!
! *** ADJUST CLOUD BASE MASS FLUX ***
!
CBMFOLD=CBMF
! *** C. Forster: adjustment of CBMF is not allowed to depend on FLEXPART timestep
DELT0=DELT/3.
DAMPS=DAMP*DELT/DELT0
CBMF=(1.-DAMPS)*CBMF+0.1*ALPHA*DTMA
CBMF=MAX(CBMF,0.0)
!
! *** If cloud base mass flux is zero, skip rest of calculation ***
!
IF(CBMF.EQ.0.0.AND.CBMFOLD.EQ.0.0)THEN
RETURN
END IF
!
! *** CALCULATE RATES OF MIXING, M(I) ***
!
M(ICB)=0.0
DO I=ICB+1,INB
K=MIN(I,INB1)
DBO=ABS(TV(K)-TVP(K))+ &
ENTP*0.02*(PHCONV_HPA(K)-PHCONV_HPA(K+1))
DBOSUM=DBOSUM+DBO
M(I)=CBMF*DBO
END DO
DO I=ICB+1,INB
M(I)=M(I)/DBOSUM
END DO
!
! *** CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING ***
! *** RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING ***
! *** FRACTION (SIJ) ***
!
DO I=ICB+1,INB
QTI=QCONV(NK)-EP(I)*CLW(I)
DO J=ICB,INB
BF2=1.+LV(J)*LV(J)*QSCONV(J)/(RV*TCONV(J)*TCONV(J)*CPD)
ANUM=H(J)-HP(I)+(CPV-CPD)*TCONV(J)*(QTI-QCONV(J))
DENOM=H(I)-HP(I)+(CPD-CPV)*(QCONV(I)-QTI)*TCONV(J)
DEI=DENOM
IF(ABS(DEI).LT.0.01)DEI=0.01
SIJ(I,J)=ANUM/DEI
SIJ(I,I)=1.0
ALTEM=SIJ(I,J)*QCONV(I)+(1.-SIJ(I,J))*QTI-QSCONV(J)
ALTEM=ALTEM/BF2
CWAT=CLW(J)*(1.-EP(J))
STEMP=SIJ(I,J)
IF((STEMP.LT.0.0.OR.STEMP.GT.1.0.OR. &
ALTEM.GT.CWAT).AND.J.GT.I)THEN
ANUM=ANUM-LV(J)*(QTI-QSCONV(J)-CWAT*BF2)
DENOM=DENOM+LV(J)*(QCONV(I)-QTI)
IF(ABS(DENOM).LT.0.01)DENOM=0.01
SIJ(I,J)=ANUM/DENOM
ALTEM=SIJ(I,J)*QCONV(I)+(1.-SIJ(I,J))*QTI-QSCONV(J)
ALTEM=ALTEM-(BF2-1.)*CWAT
END IF
IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.9)THEN
QENT(I,J)=SIJ(I,J)*QCONV(I)+(1.-SIJ(I,J))*QTI
ELIJ(I,J)=ALTEM
ELIJ(I,J)=MAX(0.0,ELIJ(I,J))
MENT(I,J)=M(I)/(1.-SIJ(I,J))
NENT(I)=NENT(I)+1
END IF
SIJ(I,J)=MAX(0.0,SIJ(I,J))
SIJ(I,J)=MIN(1.0,SIJ(I,J))
END DO
!
! *** IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS ***
! *** AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES ***
!
IF(NENT(I).EQ.0)THEN
MENT(I,I)=M(I)
QENT(I,I)=QCONV(NK)-EP(I)*CLW(I)
ELIJ(I,I)=CLW(I)
SIJ(I,I)=1.0
END IF
END DO
SIJ(INB,INB)=1.0
!
! *** NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL ***
! *** PROBABILITIES OF MIXING ***
!
DO I=ICB+1,INB
IF(NENT(I).NE.0)THEN
QP1=QCONV(NK)-EP(I)*CLW(I)
ANUM=H(I)-HP(I)-LV(I)*(QP1-QSCONV(I))
DENOM=H(I)-HP(I)+LV(I)*(QCONV(I)-QP1)
IF(ABS(DENOM).LT.0.01)DENOM=0.01
SCRIT=ANUM/DENOM
ALT=QP1-QSCONV(I)+SCRIT*(QCONV(I)-QP1)
IF(ALT.LT.0.0)SCRIT=1.0
SCRIT=MAX(SCRIT,0.0)
ASIJ=0.0
SMIN=1.0
DO J=ICB,INB
IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.9)THEN
IF(J.GT.I)THEN
SMID=MIN(SIJ(I,J),SCRIT)
SJMAX=SMID
SJMIN=SMID
IF(SMID.LT.SMIN.AND.SIJ(I,J+1).LT.SMID)THEN
SMIN=SMID
SJMAX=MIN(SIJ(I,J+1),SIJ(I,J),SCRIT)
SJMIN=MAX(SIJ(I,J-1),SIJ(I,J))
SJMIN=MIN(SJMIN,SCRIT)
END IF
ELSE
SJMAX=MAX(SIJ(I,J+1),SCRIT)
SMID=MAX(SIJ(I,J),SCRIT)
SJMIN=0.0
IF(J.GT.1)SJMIN=SIJ(I,J-1)
SJMIN=MAX(SJMIN,SCRIT)
END IF
DELP=ABS(SJMAX-SMID)
DELM=ABS(SJMIN-SMID)
ASIJ=ASIJ+(DELP+DELM)*(PHCONV_HPA(J)-PHCONV_HPA(J+1))
MENT(I,J)=MENT(I,J)*(DELP+DELM)* &
(PHCONV_HPA(J)-PHCONV_HPA(J+1))
END IF
END DO
ASIJ=MAX(1.0E-21,ASIJ)
ASIJ=1.0/ASIJ
DO J=ICB,INB
MENT(I,J)=MENT(I,J)*ASIJ
END DO
BSUM=0.0
DO J=ICB,INB
BSUM=BSUM+MENT(I,J)
END DO
IF(BSUM.LT.1.0E-18)THEN
NENT(I)=0
MENT(I,I)=M(I)
QENT(I,I)=QCONV(NK)-EP(I)*CLW(I)
ELIJ(I,I)=CLW(I)
SIJ(I,I)=1.0
END IF
END IF
END DO
!
! *** CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING ***
! *** DOWNDRAFT CALCULATION ***
!
IF(EP(INB).LT.0.0001)GOTO 405
!
! *** INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER ***
! *** AND CONDENSED WATER FLUX ***
!
JTT=2
!
! *** BEGIN DOWNDRAFT LOOP ***
!
DO I=INB,1,-1
!
! *** CALCULATE DETRAINED PRECIPITATION ***
!
WDTRAIN=G*EP(I)*M(I)*CLW(I)
IF(I.GT.1)THEN
DO J=1,I-1
AWAT=ELIJ(J,I)-(1.-EP(I))*CLW(I)
AWAT=MAX(0.0,AWAT)
WDTRAIN=WDTRAIN+G*AWAT*MENT(J,I)
END DO
END IF
!
! *** FIND RAIN WATER AND EVAPORATION USING PROVISIONAL ***
! *** ESTIMATES OF QP(I)AND QP(I-1) ***
!
!
! *** Value of terminal velocity and coefficient of evaporation for snow ***
!
COEFF=COEFFS
WT(I)=OMTSNOW
!
! *** Value of terminal velocity and coefficient of evaporation for rain ***
!
IF(TCONV(I).GT.273.0)THEN
COEFF=COEFFR
WT(I)=OMTRAIN
END IF
QSM=0.5*(QCONV(I)+QP(I+1))
AFAC=COEFF*PHCONV_HPA(I)*(QSCONV(I)-QSM)/ &
(1.0E4+2.0E3*PHCONV_HPA(I)*QSCONV(I))
AFAC=MAX(AFAC,0.0)
SIGT=SIGP(I)
SIGT=MAX(0.0,SIGT)
SIGT=MIN(1.0,SIGT)
B6=100.*(PHCONV_HPA(I)-PHCONV_HPA(I+1))*SIGT*AFAC/WT(I)
C6=(WATER(I+1)*WT(I+1)+WDTRAIN/SIGD)/WT(I)
REVAP=0.5*(-B6+SQRT(B6*B6+4.*C6))
EVAP(I)=SIGT*AFAC*REVAP
WATER(I)=REVAP*REVAP
!
! *** CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER ***
! *** HYDROSTATIC APPROXIMATION ***
!
IF(I.EQ.1)GOTO 360
DHDP=(H(I)-H(I-1))/(PCONV_HPA(I-1)-PCONV_HPA(I))
DHDP=MAX(DHDP,10.0)
MP(I)=100.*GINV*LV(I)*SIGD*EVAP(I)/DHDP
MP(I)=MAX(MP(I),0.0)
!
! *** ADD SMALL AMOUNT OF INERTIA TO DOWNDRAFT ***
!
FAC=20.0/(PHCONV_HPA(I-1)-PHCONV_HPA(I))
MP(I)=(FAC*MP(I+1)+MP(I))/(1.+FAC)
!
! *** FORCE MP TO DECREASE LINEARLY TO ZERO ***
! *** BETWEEN ABOUT 950 MB AND THE SURFACE ***
!
IF(PCONV_HPA(I).GT.(0.949*PCONV_HPA(1)))THEN
JTT=MAX(JTT,I)
MP(I)=MP(JTT)*(PCONV_HPA(1)-PCONV_HPA(I))/(PCONV_HPA(1)- &
PCONV_HPA(JTT))
END IF
360 CONTINUE
!
! *** FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT ***
!
IF(I.EQ.INB)GOTO 400
IF(I.EQ.1)THEN
QSTM=QSCONV(1)
ELSE
QSTM=QSCONV(I-1)
END IF
IF(MP(I).GT.MP(I+1))THEN
RAT=MP(I+1)/MP(I)
QP(I)=QP(I+1)*RAT+QCONV(I)*(1.0-RAT)+100.*GINV* &
SIGD*(PHCONV_HPA(I)-PHCONV_HPA(I+1))*(EVAP(I)/MP(I))
ELSE
IF(MP(I+1).GT.0.0)THEN
QP(I)=(GZ(I+1)-GZ(I)+QP(I+1)*(LV(I+1)+TCONV(I+1)*(CL-CPD))+ &
CPD*(TCONV(I+1)-TCONV(I)))/(LV(I)+TCONV(I)*(CL-CPD))
END IF
END IF
QP(I)=MIN(QP(I),QSTM)
QP(I)=MAX(QP(I),0.0)
400 CONTINUE
END DO
!
! *** CALCULATE SURFACE PRECIPITATION IN MM/DAY ***
!
PRECIP=PRECIP+WT(1)*SIGD*WATER(1)*3600.*24000./(ROWL*G)
!
405 CONTINUE
!
! *** CALCULATE DOWNDRAFT VELOCITY SCALE AND SURFACE TEMPERATURE AND ***
! *** WATER VAPOR FLUCTUATIONS ***
!
WD=BETA*ABS(MP(ICB))*0.01*RD*TCONV(ICB)/(SIGD*PCONV_HPA(ICB))
QPRIME=0.5*(QP(1)-QCONV(1))
TPRIME=LV0*QPRIME/CPD
!
! *** CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE ***
! *** AND MIXING RATIO ***
!
DPINV=0.01/(PHCONV_HPA(1)-PHCONV_HPA(2))
AM=0.0
IF(NK.EQ.1)THEN
DO K=2,INB
AM=AM+M(K)
END DO
END IF
! save saturated upward mass flux for first level
FUP(1)=AM
IF((2.*G*DPINV*AM).GE.DELTI)IFLAG=4
FT(1)=FT(1)+G*DPINV*AM*(TCONV(2)-TCONV(1)+(GZ(2)-GZ(1))/CPN(1))
FT(1)=FT(1)-LVCP(1)*SIGD*EVAP(1)
FT(1)=FT(1)+SIGD*WT(2)*(CL-CPD)*WATER(2)*(TCONV(2)- &
TCONV(1))*DPINV/CPN(1)
FQ(1)=FQ(1)+G*MP(2)*(QP(2)-QCONV(1))* &
DPINV+SIGD*EVAP(1)
FQ(1)=FQ(1)+G*AM*(QCONV(2)-QCONV(1))*DPINV
DO J=2,INB
FQ(1)=FQ(1)+G*DPINV*MENT(J,1)*(QENT(J,1)-QCONV(1))
END DO
!
! *** CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO ***
! *** AT LEVELS ABOVE THE LOWEST LEVEL ***
!
! *** FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES ***
! *** THROUGH EACH LEVEL ***
!
DO I=2,INB
DPINV=0.01/(PHCONV_HPA(I)-PHCONV_HPA(I+1))
CPINV=1.0/CPN(I)
AMP1=0.0
AD=0.0
IF(I.GE.NK)THEN
DO K=I+1,INB+1
AMP1=AMP1+M(K)
END DO
END IF
DO K=1,I
DO J=I+1,INB+1
AMP1=AMP1+MENT(K,J)
END DO
END DO
! save saturated upward mass flux
FUP(I)=AMP1
IF((2.*G*DPINV*AMP1).GE.DELTI)IFLAG=4
DO K=1,I-1
DO J=I,INB
AD=AD+MENT(J,K)
END DO
END DO
! save saturated downward mass flux
FDOWN(I)=AD
FT(I)=FT(I)+G*DPINV*(AMP1*(TCONV(I+1)-TCONV(I)+(GZ(I+1)-GZ(I))* &
CPINV)-AD*(TCONV(I)-TCONV(I-1)+(GZ(I)-GZ(I-1))*CPINV)) &
-SIGD*LVCP(I)*EVAP(I)
FT(I)=FT(I)+G*DPINV*MENT(I,I)*(HP(I)-H(I)+ &
TCONV(I)*(CPV-CPD)*(QCONV(I)-QENT(I,I)))*CPINV
FT(I)=FT(I)+SIGD*WT(I+1)*(CL-CPD)*WATER(I+1)* &
(TCONV(I+1)-TCONV(I))*DPINV*CPINV
FQ(I)=FQ(I)+G*DPINV*(AMP1*(QCONV(I+1)-QCONV(I))- &
AD*(QCONV(I)-QCONV(I-1)))
DO K=1,I-1
AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I)
AWAT=MAX(AWAT,0.0)
FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-AWAT-QCONV(I))
END DO
DO K=I,INB
FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-QCONV(I))
END DO
FQ(I)=FQ(I)+SIGD*EVAP(I)+G*(MP(I+1)* &
(QP(I+1)-QCONV(I))-MP(I)*(QP(I)-QCONV(I-1)))*DPINV
END DO
!
! *** Adjust tendencies at top of convection layer to reflect ***
! *** actual position of the level zero CAPE ***
!
FQOLD=FQ(INB)
FQ(INB)=FQ(INB)*(1.-FRAC)
FQ(INB-1)=FQ(INB-1)+FRAC*FQOLD*((PHCONV_HPA(INB)- &
PHCONV_HPA(INB+1))/ &
(PHCONV_HPA(INB-1)-PHCONV_HPA(INB)))*LV(INB)/LV(INB-1)
FTOLD=FT(INB)
FT(INB)=FT(INB)*(1.-FRAC)
FT(INB-1)=FT(INB-1)+FRAC*FTOLD*((PHCONV_HPA(INB)- &
PHCONV_HPA(INB+1))/ &
(PHCONV_HPA(INB-1)-PHCONV_HPA(INB)))*CPN(INB)/CPN(INB-1)
!
! *** Very slightly adjust tendencies to force exact ***
! *** enthalpy, momentum and tracer conservation ***
!
ENTS=0.0
DO I=1,INB
ENTS=ENTS+(CPN(I)*FT(I)+LV(I)*FQ(I))* &
(PHCONV_HPA(I)-PHCONV_HPA(I+1))
END DO
ENTS=ENTS/(PHCONV_HPA(1)-PHCONV_HPA(INB+1))
DO I=1,INB
FT(I)=FT(I)-ENTS/CPN(I)
END DO
! ************************************************
! **** DETERMINE MASS DISPLACEMENT MATRIX
! ***** AND COMPENSATING SUBSIDENCE
! ************************************************
! mass displacement matrix due to saturated up-and downdrafts
! inside the cloud and determine compensating subsidence
! FUP(I) (saturated updrafts), FDOWN(I) (saturated downdrafts) are assumed to be
! balanced by compensating subsidence (SUB(I))
! FDOWN(I) and SUB(I) defined positive downwards
! NCONVTOP IS THE TOP LEVEL AT WHICH CONVECTIVE MASS FLUXES ARE DIAGNOSED
! EPSILON IS A SMALL NUMBER
SUB(1)=0.
NCONVTOP=1
do i=1,INB+1
do j=1,INB+1
if (j.eq.NK) then
FMASS(j,i)=FMASS(j,i)+M(i)
endif
FMASS(j,i)=FMASS(j,i)+MENT(j,i)
IF (FMASS(J,I).GT.EPSILON) NCONVTOP=MAX(NCONVTOP,I,J)
end do
if (i.gt.1) then
SUB(i)=FUP(i-1)-FDOWN(i)
endif
end do
NCONVTOP=NCONVTOP+1
RETURN
!
END SUBROUTINE CONVECT
!
! ---------------------------------------------------------------------------
!
SUBROUTINE TLIFT(GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK)
!
!-cv
use par_mod
use conv_mod
implicit none
!-cv
!====>Begin Module TLIFT File convect.f Undeclared variables
!
!Argument variables
!
integer :: icb, kk, nd, nk, nl
!
!Local variables
!
integer :: i, j, nsb, nst
!
real :: ah0, ahg, alv, cpinv, cpp, denom
real :: es, qg, rg, s, tc, tg
!
!====>End Module TLIFT File convect.f
REAL :: GZ(ND),TPK(ND),CLW(ND)
REAL :: TVP(ND)
!
! *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS ***
!
REAL,PARAMETER :: CPD=1005.7
REAL,PARAMETER :: CPV=1870.0
REAL,PARAMETER :: CL=2500.0
REAL,PARAMETER :: RV=461.5
REAL,PARAMETER :: RD=287.04
REAL,PARAMETER :: LV0=2.501E6
!
REAL,PARAMETER :: CPVMCL=CL-CPV
REAL,PARAMETER :: EPS0=RD/RV
REAL,PARAMETER :: EPSI=1./EPS0
!
! *** CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY ***
!
AH0=(CPD*(1.-QCONV(NK))+CL*QCONV(NK))*TCONV(NK)+QCONV(NK)* &
(LV0-CPVMCL*( &
TCONV(NK)-273.15))+GZ(NK)
CPP=CPD*(1.-QCONV(NK))+QCONV(NK)*CPV
CPINV=1./CPP
!
IF(KK.EQ.1)THEN
!
! *** CALCULATE LIFTED PARCEL QUANTITIES BELOW CLOUD BASE ***
!
DO I=1,ICB-1
CLW(I)=0.0
END DO
DO I=NK,ICB-1
TPK(I)=TCONV(NK)-(GZ(I)-GZ(NK))*CPINV
TVP(I)=TPK(I)*(1.+QCONV(NK)*EPSI)
END DO
END IF
!
! *** FIND LIFTED PARCEL QUANTITIES ABOVE CLOUD BASE ***
!
NST=ICB
NSB=ICB
IF(KK.EQ.2)THEN
NST=NL
NSB=ICB+1
END IF
DO I=NSB,NST
TG=TCONV(I)
QG=QSCONV(I)
ALV=LV0-CPVMCL*(TCONV(I)-273.15)
DO J=1,2
S=CPD+ALV*ALV*QG/(RV*TCONV(I)*TCONV(I))
S=1./S
AHG=CPD*TG+(CL-CPD)*QCONV(NK)*TCONV(I)+ALV*QG+GZ(I)
TG=TG+S*(AH0-AHG)
TG=MAX(TG,35.0)
TC=TG-273.15
DENOM=243.5+TC
IF(TC.GE.0.0)THEN
ES=6.112*EXP(17.67*TC/DENOM)
ELSE
ES=EXP(23.33086-6111.72784/TG+0.15215*LOG(TG))
END IF
QG=EPS0*ES/(PCONV_HPA(I)-ES*(1.-EPS0))
END DO
ALV=LV0-CPVMCL*(TCONV(I)-273.15)
TPK(I)=(AH0-(CL-CPD)*QCONV(NK)*TCONV(I)-GZ(I)-ALV*QG)/CPD
CLW(I)=QCONV(NK)-QG
CLW(I)=MAX(0.0,CLW(I))
RG=QG/(1.-QCONV(NK))
TVP(I)=TPK(I)*(1.+RG*EPSI)
END DO
RETURN
END SUBROUTINE TLIFT
|
