forked from wrf-model/WRF
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmodule_bl_camuwpbl_driver.F
executable file
·1097 lines (910 loc) · 63.1 KB
/
module_bl_camuwpbl_driver.F
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
#define MODAL_AERO
MODULE module_bl_camuwpbl_driver
!Note: Comments starting with "!!" are directly taken from CAM's interface routine for the UW PBL (vertical_diffusion.F90)
!This modules is based on vertical_diffusion.F90 of CAM5
!List of possible modifications in the future:
!_____________________________________________
!1. Constituents(per Kg of DRY air or "dry" constituents) are not diffused currently. The call which
! makes DRY constituents diffusion decision flag 'true' is commented out in camuwpblinit subroutine.
! Therefore the flag is ALWAYS false for DRY constituent diffusion.
!
! REASON: DRY constituents are Aerosols and other species (excepts for cloud mass mixing ratios and number concentrations).
! ------ In WRF, DRY constituents appear only when WRF_CHEM simulations are run. DRY constituents are vertically mixed(diffused)
! in dry_dep_driver of WRF_CHEM. Therefore, DRY constituents are not diffused in this PBL subroutine
!
!2. The liquid number concentrations is NOT diffused(or mixed) in PBL to mimic CAM5, which diffuses liquid number
! concentrations in dropmixnuc.
!3. Surface fluxes for ALL the constituents are ZERO except for the water vapours (1st constituent)
!4. Molecular diffusion is turned off for now
!5. Mountain stresses are not computed for now (all calls to Mountain stresses currently may hold undefined variables)
!6. *Formulas used for computing surface stresses is based on the formula given in CAM5's BareGround.F90.
! This formula may not work well for ocean surfaces (possible future modification)
!_____________________________________________
!!----------------------------------------------------------------------------------------------------- !
!! Module to compute vertical diffusion of momentum, moisture, trace constituents !
!! and static energy. Separate modules compute !
!! 1. stresses associated with turbulent flow over orography !
!! ( turbulent mountain stress ) !
!! 2. eddy diffusivities, including nonlocal tranport terms !
!! 3. molecular diffusivities !
!! 4. coming soon... gravity wave drag !
!! Lastly, a implicit diffusion solver is called, and tendencies retrieved by !
!! differencing the diffused and initial states. !
!! !
!! Calling sequence: !
!! !
!! vertical_diffusion_init Initializes vertical diffustion constants and modules !
!! init_molec_diff Initializes molecular diffusivity module !
!! init_eddy_diff Initializes eddy diffusivity module (includes PBL) !
!! init_tms Initializes turbulent mountain stress module !
!! init_vdiff Initializes diffusion solver module !
!! vertical_diffusion_ts_init Time step initialization (only used for upper boundary condition) !
!! vertical_diffusion_tend Computes vertical diffusion tendencies !
!! compute_tms Computes turbulent mountain stresses !
!! compute_eddy_diff Computes eddy diffusivities and countergradient terms !
!! compute_vdiff Solves vertical diffusion equations, including molecular diffusivities !
!! !
!!---------------------------Code history-------------------------------------------------------------- !
!! J. Rosinski : Jun. 1992 !
!! J. McCaa : Sep. 2004 !
!! S. Park : Aug. 2006, Dec. 2008. Jan. 2010 !
! B. Singh : Nov. 2010 (ported to WRF by balwinder.singh@pnl.gov)
!!----------------------------------------------------------------------------------------------------- !
use shr_kind_mod, only : r8 => shr_kind_r8
use module_cam_support, only : pcols, pver, pverp, endrun, iulog,fieldname_len,pcnst =>pcnst_runtime
use constituents, only : qmin
use diffusion_solver, only : vdiff_selector
use physconst, only : &
cpair , & ! Specific heat of dry air
gravit , & ! Acceleration due to gravity
rair , & ! Gas constant for dry air
zvir , & ! rh2o/rair - 1
latvap , & ! Latent heat of vaporization
latice , & ! Latent heat of fusion
karman , & ! von Karman constant
mwdry , & ! Molecular weight of dry air
avogad , & ! Avogadro's number
boltz ! Boltzman's constant
implicit none
private
save
!! ----------------- !
!! Public interfaces !
!! ----------------- !
public camuwpblinit ! Initialization
public camuwpbl ! Driver for the PBL scheme
public vd_register ! Init routine for constituents
!! ------------ !
!! Private data !
!! ------------ !
character(len=16) :: eddy_scheme !! Default set in phys_control.F90, use namelist to change
!! 'HB' = Holtslag and Boville (default)
!! 'HBR' = Holtslag and Boville and Rash
!! 'diag_TKE' = Bretherton and Park ( UW Moist Turbulence Scheme )
integer, parameter :: nturb = 5 !! Number of iterations for solution ( when 'diag_TKE' scheme is selected )
logical, parameter :: wstarent = .true. !! Use wstar (.true.) or TKE (.false.) entrainment closure ( when 'diag_TKE' scheme is selected )
logical :: do_pseudocon_diff = .false. !! If .true., do pseudo-conservative variables diffusion
character(len=16) :: shallow_scheme !! For checking compatibility between eddy diffusion and shallow convection schemes
!! 'Hack' = Hack Shallow Convection Scheme
!! 'UW' = Park and Bretherton ( UW Shallow Convection Scheme )
character(len=16) :: microp_scheme !! Microphysics scheme
logical :: do_molec_diff = .false. !! Switch for molecular diffusion
logical :: do_tms !! Switch for turbulent mountain stress
real(r8) :: tms_orocnst !! Converts from standard deviation to height
real(r8) :: tms_z0fac ! Converts from standard deviation to height
type(vdiff_selector) :: fieldlist_wet !! Logical switches for moist mixing ratio diffusion
type(vdiff_selector) :: fieldlist_dry !! Logical switches for dry mixing ratio diffusion
integer :: ntop !! Top interface level to which vertical diffusion is applied ( = 1 ).
integer :: nbot !! Bottom interface level to which vertical diffusion is applied ( = pver ).
integer :: tke_idx, kvh_idx, kvm_idx !! TKE and eddy diffusivity indices for fields in the physics buffer
integer :: turbtype_idx, smaw_idx !! Turbulence type and instability functions
integer :: tauresx_idx, tauresy_idx !! Redisual stress for implicit surface stress
integer :: ixcldice, ixcldliq !! Constituent indices for cloud liquid and ice water
integer :: ixnumice, ixnumliq
#ifdef MODAL_AERO
integer :: ixndrop
#endif
integer :: wgustd_index
CONTAINS
subroutine camuwpbl(dt,u_phy,v_phy,th_phy,rho,qv_curr,hfx,qfx,ustar,p8w &
,p_phy,z,t_phy,qc_curr,qi_curr,z_at_w,cldfra_old_mp,cldfra,ht &
,rthratenlw,exner,is_CAMMGMP_used &
,itimestep,qnc_curr,qni_curr,wsedl3d &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
!Intent-inout
,tauresx2d,tauresy2d &
,rublten,rvblten,rthblten,rqiblten,rqniblten,rqvblten,rqcblten &
,kvm3d,kvh3d &
!Intent-out
,tpert2d,qpert2d,wpert2d,smaw3d,turbtype3d &
,tke_pbl,pblh2d,kpbl2d )
!!---------------------------------------------------- !
!! This is an interface routine for vertical diffusion !
! This subroutine is called by module_pbl_driver and !
! it calls: wrf_error_fatal,compute_tms, !
! compute_eddy_diff,aqsat,compute_vdiff and !
! positive_moisture.
!!---------------------------------------------------- !
use module_cam_support, only : pcols
use trb_mtn_stress, only : compute_tms
use eddy_diff, only : compute_eddy_diff
use wv_saturation, only : fqsatd, aqsat
use molec_diff, only : compute_molec_diff, vd_lu_qdecomp
use constituents, only : qmincg, qmin
use diffusion_solver !!, only : compute_vdiff, any, operator(.not.)
#ifdef MODAL_AERO
use modal_aero_data
#endif
implicit none
!------------------------------------------------------------------------!
! Input !
!------------------------------------------------------------------------!
logical, intent(in) :: is_CAMMGMP_used
integer, intent(in) :: ids,ide, jds,jde, kds,kde
integer, intent(in) :: ims,ime, jms,jme, kms,kme
integer, intent(in) :: its,ite, jts,jte, kts,kte
integer, intent(in) :: itimestep
real, intent(in) :: dt ! Time step (s)
real, dimension( ims:ime,jms:jme ), intent(in) :: hfx !Sensible heat flux at surface (w/m2)
real, dimension( ims:ime,jms:jme ), intent(in) :: qfx !Moisture flux at surface (kg m-2 s-1)
real, dimension( ims:ime,jms:jme ), intent(in) :: ustar !Friction velocity (m/s)
real, dimension( ims:ime,jms:jme ), intent(in) :: ht !Terrain height (m)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: v_phy !Y-component of wind (m/s)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: u_phy !X-component of wind (m/s)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: th_phy !Potential temperature (K)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: rho !Air density (kg/m3)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qv_curr !Water vapor mixing ratio - Moisture (kg/kg)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qc_curr !Cloud water mixing ratio - Cloud liq (kg/kg)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qi_curr !Ice mixing ratio - Cloud ice (kg/kg)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qnc_curr !Liq # mixing ratio - Cloud liq # (#/kg)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qni_curr !Ice # mixing ratio - Cloud ice # (#/kg)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p_phy !Pressure at mid-level (Pa)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p8w !Pressure at level interface (Pa)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z !Height above sea level at mid-level (m)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: t_phy !Temperature (K)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z_at_w !Height above sea level at layer interfaces (m)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cldfra_old_mp !Cloud fraction [unitless]
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cldfra !Cloud fraction [unitless]
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: exner !Dimensionless pressure [unitless]
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: rthratenlw !Tendency for LW ( K/s)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: wsedl3d !Sedimentation velocity of stratiform liquid cloud droplet (m/s)
!------------------------------------------------------------------------!
! Output !
!------------------------------------------------------------------------!
real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rublten !Tendency for u_phy (Pa m s-2)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rvblten !Tendency for v_phy (Pa m s-2)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rthblten !Tendency for th_phy (Pa K s-1)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqvblten !Tendency for qv_curr (Pa kg kg-1 s-1)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqcblten !Tendency for qc_curr (Pa kg kg-1 s-1)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqiblten !Tendency for qi_curr (Pa kg kg-1 s-1)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqniblten !Tendency for qni_curr(Pa # kg-1 s-1)
integer, dimension( ims:ime,jms:jme ), intent(out) :: kpbl2d !Layer index containing PBL top within or at the base interface
real, dimension( ims:ime,jms:jme ), intent(out) :: pblh2d !Planetary boundary layer height (m)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: tke_pbl !Turbulence kinetic energy at midpoints (m^2/s^2)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: turbtype3d !Turbulent interface types [ no unit ]
real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: smaw3d !Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units]
!---------------------------------------------------------------------------!
! Local, carried on from one timestep to the other (defined in registry)!
!---------------------------------------------------------------------------!
real, dimension( ims:ime, jms:jme ) , intent(inout ):: tauresx2d,tauresy2d !X AND Y-COMP OF RESIDUAL STRESSES(m^2/s^2)
real, dimension( ims:ime, jms:jme ) , intent(out) :: tpert2d ! Convective temperature excess (K)
real, dimension( ims:ime, jms:jme ) , intent(out) :: qpert2d ! Convective humidity excess (kg/kg)
real, dimension( ims:ime, jms:jme ) , intent(out) :: wpert2d ! Turbulent velocity excess (m/s)
real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: kvm3d,kvh3d !Eddy diffusivity for momentum and heat(m^2/s)
!---------------------------------------------------------------------------!
! Local !
!---------------------------------------------------------------------------!
character(128) :: errstring ! Error status for compute_vdiff
logical :: kvinit ! Tell compute_eddy_diff/ caleddy to initialize kvh, kvm (uses kvf)
logical :: is_first_step ! Flag to know if this a first time step
integer :: i,j,k,itsp1,itile_len,ktep1,kflip,ncol,ips,icnst,m,kp1
integer :: lchnk ! Chunk identifier
real(r8) :: tauFac, uMean, dp, multFrc
real(r8) :: ztodt ! 2*delta-t (s)
real(r8) :: rztodt ! 1./ztodt (1/s)
real(r8),parameter :: invalidVal = -999888777.0_r8
real(r8) :: topflx( pcols) ! Molecular heat flux at top interface
real(r8) :: wpert( pcols) ! Turbulent velocity excess (m/s)
real(r8) :: tauresx( pcols) ! [Residual stress to be added in vdiff to correct...
real(r8) :: tauresy( pcols) ! for turb stress mismatch between sfc and atm accumulated.] (N/m2)
real(r8) :: ipbl( pcols) ! If 1, PBL is CL, while if 0, PBL is STL.
real(r8) :: kpblh( pcols) ! Layer index containing PBL top within or at the base interface
real(r8) :: wstarPBL(pcols) ! Convective velocity within PBL (m/s)
real(r8) :: sgh( pcols) ! Standard deviation of orography (m)
real(r8) :: landfrac(pcols) ! Land fraction [unitless]
real(r8) :: taux( pcols) ! x surface stress (N/m2)
real(r8) :: tauy( pcols) ! y surface stress (N/m2)
real(r8) :: tautotx( pcols) ! U component of total surface stress (N/m2)
real(r8) :: tautoty( pcols) ! V component of total surface stress (N/m2)
real(r8) :: ksrftms( pcols) ! Turbulent mountain stress surface drag coefficient (kg/s/m2)
real(r8) :: tautmsx( pcols) ! U component of turbulent mountain stress (N/m2)
real(r8) :: tautmsy( pcols) ! V component of turbulent mountain stress (N/m2)
real(r8) :: ustar8( pcols) ! Surface friction velocity (m/s)
real(r8) :: pblh( pcols) ! Planetary boundary layer height (m)
real(r8) :: tpert( pcols) ! Convective temperature excess (K)
real(r8) :: qpert( pcols) ! Convective humidity excess (kg/kg)
real(r8) :: shflx( pcols) ! Surface sensible heat flux (w/m2)
real(r8) :: phis( pcols) ! Geopotential at terrain height (m2/s2)
real(r8) :: cldn8( pcols,kte) ! New stratus fraction (fraction)
real(r8) :: qrl8( pcols,kte) ! LW radiative cooling rate(W/kg*Pa)
real(r8) :: wsedl8( pcols,kte) ! Sedimentation velocity of stratiform liquid cloud droplet (m/s)
real(r8) :: dtk( pcols,kte) ! T tendency from KE dissipation
real(r8) :: qt( pcols,kte) !
real(r8) :: sl_prePBL( pcols,kte) !
real(r8) :: qt_prePBL( pcols,kte) !
real(r8) :: slten( pcols,kte) !
real(r8) :: qtten( pcols,kte) !
real(r8) :: sl( pcols,kte) !
real(r8) :: ftem( pcols,kte) ! Saturation vapor pressure before PBL
real(r8) :: ftem_prePBL(pcols,kte) ! Saturation vapor pressure before PBL
real(r8) :: ftem_aftPBL(pcols,kte) ! Saturation vapor pressure after PBL
real(r8) :: tem2( pcols,kte) ! Saturation specific humidity and RH
real(r8) :: t_aftPBL( pcols,kte) ! Temperature after PBL diffusion
real(r8) :: tten( pcols,kte) ! Temperature tendency by PBL diffusion
real(r8) :: rhten( pcols,kte) ! RH tendency by PBL diffusion
real(r8) :: qv_aft_PBL( pcols,kte) ! qv after PBL diffusion
real(r8) :: ql_aft_PBL( pcols,kte) ! ql after PBL diffusion
real(r8) :: qi_aft_PBL( pcols,kte) ! qi after PBL diffusion
real(r8) :: s_aft_PBL( pcols,kte) ! s after PBL diffusion
real(r8) :: u_aft_PBL( pcols,kte) ! u after PBL diffusion
real(r8) :: v_aft_PBL( pcols,kte) ! v after PBL diffusion
real(r8) :: qv_pro( pcols,kte) !
real(r8) :: ql_pro( pcols,kte) !
real(r8) :: qi_pro( pcols,kte) !
real(r8) :: s_pro( pcols,kte) !
real(r8) :: t_pro( pcols,kte) !
real(r8) :: u8( pcols,kte) ! x component of velocity in CAM's data structure (m/s)
real(r8) :: v8( pcols,kte) ! y component of velocity in CAM's data structure (m/s)
real(r8) :: t8( pcols,kte) !
real(r8) :: pmid8( pcols,kte) ! Pressure at the midpoints in CAM's data structure (Pa)
real(r8) :: pmiddry8( pcols,kte) ! Dry Pressure at the midpoints in CAM's data structure (Pa)
real(r8) :: zm8( pcols,kte) ! Height at the midpoints in CAM's data structure (m)
real(r8) :: exner8( pcols,kte) ! exner function in CAM's data structure
real(r8) :: s8( pcols,kte) ! Dry static energy (m2/s2)
real(r8) :: rpdel8( pcols,kte) ! Inverse of pressure difference (1/Pa)
real(r8) :: pdel8( pcols,kte) ! Pressure difference (Pa)
real(r8) :: rpdeldry8( pcols,kte) ! Inverse of dry pressure difference (1/Pa)
real(r8) :: pdeldry8( pcols,kte) ! dry pressure difference (1/Pa)
REAL(r8) :: stnd( pcols,kte) ! Heating rate (dry static energy tendency, W/kg)
real(r8) :: tke8( pcols,kte+1) ! Turbulent kinetic energy [ m2/s2 ]
real(r8) :: turbtype( pcols,kte+1) ! Turbulent interface types [ no unit ]
real(r8) :: smaw( pcols,kte+1) ! Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units]
! This is 1 when neutral condition (Ri=0), 4.964 for maximum unstable case, and 0 when Ri > Ricrit=0.19.
real(r8) :: cgs( pcols,kte+1) ! Counter-gradient star [ cg/flux ]
real(r8) :: cgh( pcols,kte+1) ! Counter-gradient term for heat [ J/kg/m ]
real(r8) :: kvh( pcols,kte+1) ! Eddy diffusivity for heat [ m2/s ]
real(r8) :: kvm( pcols,kte+1) ! Eddy diffusivity for momentum [ m2/s ]
real(r8) :: kvq( pcols,kte+1) ! Eddy diffusivity for constituents [ m2/s ]
real(r8) :: kvh_in( pcols,kte+1) ! kvh from previous timestep [ m2/s ]
real(r8) :: kvm_in( pcols,kte+1) ! kvm from previous timestep [ m2/s ]
real(r8) :: bprod( pcols,kte+1) ! Buoyancy production of tke [ m2/s3 ]
real(r8) :: sprod( pcols,kte+1) ! Shear production of tke [ m2/s3 ]
real(r8) :: sfi( pcols,kte+1) ! Saturation fraction at interfaces [ fraction ]
real(r8) :: pint8( pcols,kte+1) ! Pressure at the interfaces in CAM's data structure (Pa)
real(r8) :: pintdry8( pcols,kte+1) ! Dry pressure at the interfaces in CAM's data structure (Pa)
real(r8) :: zi8( pcols,kte+1) ! Height at the interfacesin CAM's data structure (m)
real(r8) :: cloud( pcols,kte,pcnst) ! Holder for cloud water and ice (q in cam)
real(r8) :: cloudtnd( pcols,kte,pcnst) ! Holder for cloud tendencies (ptend_loc%q in cam)
real(r8) :: wind_tends(pcols,kte,2) ! Wind component tendencies (m/s2)
real(r8) :: cflx(pcols,pcnst) ! Surface constituent flux [ kg/m2/s ]
#ifdef MODAL_AERO
real(r8) :: tmp1(pcols) ! Temporary storage
integer :: l, lspec
#endif
!! ----------------------- !
!! Main Computation Begins !
!! ----------------------- !
do icnst = 1 , pcnst
!Setting curface fluxes for all constituents to be zero. Later, cflx(:,1) is populated by the water vapour surface flux
cflx(:pcols,icnst) = 0.0_r8
enddo
is_first_step = .false.
if(itimestep == 1) then
is_first_step = .true.
endif
ncol = pcols
ztodt = DT
rztodt = 1.0_r8 / ztodt
!Few definitions in this subroutine require that ncol==1
if(ncol .NE. 1) then
call wrf_error_fatal('Number of CAM Columns (NCOL) in CAMUWPBL scheme must be 1')
endif
!Initialize errstring
errstring = ''
!-------------------------------------------------------------------------------------!
!Map CAM variables to the corresponding WRF variables !
!Loop over the points in the tile and treat them each as a CAM Chunk !
!-------------------------------------------------------------------------------------!
itsp1 = its - 1
itile_len = ite - itsp1
do j = jts , jte
do i = its , ite
lchnk = (j - jts) * itile_len + (i - itsp1) !1-D index location from a 2-D tile
phis(1) = ht(i,j) * gravit !Used for computing dry static energy
!Flip vertically quantities computed at the mid points
ktep1 = kte + 1
do k = kts,kte
kflip = ktep1 - k
!Loop is used as vector assignment may take more computational time
do icnst = 1 , pcnst
!Setting cloud and cloudtnd to values (obtained from module_cam_support) which should produce errors if used
!1st,2nd,3rd and 5th constituents are used (see the assignments below) and diffused presently in this scheme
cloud(1,kflip,icnst) = invalidVal
cloudtnd(1,kflip,icnst) = invalidVal
enddo
u8( 1,kflip) = u_phy(i,k,j) ! X-component of velocity at the mid-points (m/s) [state%u in CAM]
v8( 1,kflip) = v_phy(i,k,j) ! Y-component of velocity at the mid-points (m/s) [state%v in CAM]
pmid8( 1,kflip) = p_phy(i,k,j) ! Pressure at the mid-points (Pa) [state%pmid in CAM]
dp = p8w(i,k,j) - p8w(i,k+1,j) !Change in pressure (Pa)
pdel8( 1,kflip) = dp
rpdel8( 1,kflip) = 1.0_r8/dp ! Reciprocal of pressure difference (1/Pa) [state%rpdel in CAM]
zm8( 1,kflip) = z(i,k,j)-ht(i,j) ! Height above the ground at the midpoints (m) [state%zm in CAM]
t8( 1,kflip) = t_phy(i,k,j) ! Temprature at the mid points (K) [state%t in CAM]
s8( 1,kflip) = cpair *t8(1,kflip) + gravit*zm8(1,kflip) + phis(1) ! Dry static energy (m2/s2) [state%s in CAM]-Formula tested in vertical_diffusion.F90 of CAM
qrl8( 1,kflip) = rthratenlw(i,k,j) * exner(i,k,j) * cpair * dp ! Long Wave heating rate (W/kg*Pa)- Formula obtained from definition of qrlw(pcols,pver) in eddy_diff.F90 in CAM
wsedl8( 1,kflip) = wsedl3d(i,k,j) ! Sedimentation velocity of stratiform liquid cloud droplet (m/s)
!Following three formulas are obtained from ported CAM's ZM cumulus scheme
!Values of 0 cause a crash in entropy
multFrc = 1._r8/(1._r8 + qv_curr(i,k,j))
cloud( 1,kflip,1) = max( qv_curr(i,k,j) * multFrc, 1.0e-30_r8 ) !Specific humidity [state%q(:,:,1) in CAM]
cloud( 1,kflip,ixcldliq) = qc_curr(i,k,j) * multFrc !Convert to moist mix ratio-cloud liquid [state%q(:,:,2) in CAM]
cloud( 1,kflip,ixcldice) = qi_curr(i,k,j) * multFrc !cloud ice [state%q(:,:,3) in CAM]
cloud( 1,kflip,ixnumliq) = qnc_curr(i,k,j) * multFrc
cloud( 1,kflip,ixnumice) = qni_curr(i,k,j) * multFrc
exner8(1,kflip) = exner(i,k,j) !Exner function (no units)
if(is_CAMMGMP_used) then
cldn8( 1,kflip) = cldfra_old_mp(i,k,j) !Cloud fraction (no unit)
else
cldn8( 1,kflip) = cldfra(i,k,j) !Cloud fraction (no unit)
endif
!Following formula is obtained from physics_types.F90 of CAM (CESM101)
pdeldry8(1,kflip) = pdel8(1,kflip)*(1._r8-cloud(1,kflip,1))
rpdeldry8(1,kflip) = 1._r8/pdeldry8(1,kflip)
enddo
do k = kts,kte+1
kflip = kte - k + 2
pint8( 1,kflip) = p8w( i,k,j) ! Pressure at interfaces [state%pint in CAM]
zi8( 1,kflip) = z_at_w(i,k,j) -ht(i,j) ! Height at interfaces [state%zi in CAM]
!Initialize Variables to zero, these are outputs from the "compute_eddy_diff"
kvq(1,kflip) = 0.0_r8 ! Eddy diffusivity for constituents (m2/s)
cgh(1,kflip) = 0.0_r8 ! Counter-gradient term for heat
cgs(1,kflip) = 0.0_r8 ! Counter-gradient star (cg/flux)
if( is_first_step ) then
kvm3d(i,k,j) = 0.0_r8 ! Eddy diffusivity for heat (m2/s)
kvh3d(i,k,j) = 0.0_r8 ! Eddy diffusivity for momentum (m2/s)
endif
kvh(1,kflip) = kvh3d(i,k,j)
kvm(1,kflip) = kvm3d(i,k,j)
end do
!Compute pintdry8 and pmiddry8
!Following formula is obtained from physics_types.F90 of CAM (CESM101)
pintdry8(1,1) = pint8(1,1)
do k = 1, pver
pintdry8(1,k+1) = pintdry8(1,k) + pdeldry8(1,k)
pmiddry8(1,k) = (pintdry8(1,k+1)+ pintdry8(1,k))*0.5_r8
end do
shflx( ncol) = hfx(i,j) ! Surface sensible heat flux (w/m2)
!SGH and LANDFRAC are inputs for the compute_tms subroutine. Presently set to zero as do_tms is always false for now
sgh( ncol) = 0.0_r8 ! Standard deviation of orography (m)
landfrac(ncol) = 0.0_r8 ! Fraction (unitless)
uMean = sqrt( u_phy(i,kts,j) * u_phy(i,kts,j) + v_phy(i,kts,j) * v_phy(i,kts,j) ) ! Mean velocity
tauFac = rho(i,kts,j) * ustar(i,j) *ustar(i,j)/uMean
taux(ncol) = -tauFac * u_phy(i,kts,j) ! x surface stress (N/m2) [Formulation obtained from CAM's BareGround.F90]
tauy(ncol) = -tauFac * v_phy(i,kts,j) ! y surface stress (N/m2)
!! Retrieve 'tauresx, tauresy' from from the last timestep
if( is_first_step ) then
tauresx2d(i,j) = 0._r8
tauresy2d(i,j) = 0._r8
endif
tauresx(:ncol) = tauresx2d(i,j)
tauresy(:ncol) = tauresy2d(i,j)
!! All variables are modified by vertical diffusion
!!------------------------------------------!
!! Computation of turbulent mountain stress !
!!------------------------------------------!
!! Consistent with the computation of 'normal' drag coefficient, we are using
!! the raw input (u,v) to compute 'ksrftms', not the provisionally-marched 'u,v'
!! within the iteration loop of the PBL scheme.
if( do_tms ) then
call compute_tms( pcols , pver , ncol , &
u8 , v8 , t8 , pmid8 , &
exner8 , zm8 , sgh , ksrftms , &
tautmsx , tautmsy , landfrac )
!! Here, both 'taux, tautmsx' are explicit surface stresses.
!! Note that this 'tautotx, tautoty' are different from the total stress
!! that has been actually added into the atmosphere. This is because both
!! taux and tautmsx are fully implicitly treated within compute_vdiff.
!! However, 'tautotx, tautoty' are not used in the actual numerical
!! computation in this module.
tautotx(:ncol) = taux(:ncol) + tautmsx(:ncol)
tautoty(:ncol) = tauy(:ncol) + tautmsy(:ncol)
else
ksrftms(:ncol) = 0.0_r8
tautotx(:ncol) = taux(:ncol)
tautoty(:ncol) = tauy(:ncol)
endif
!-------------------------------------------------------------------------------------!
!We are currenly using just water vapour flux at the surface, rest are set to zero
!-------------------------------------------------------------------------------------!
cflx(:pcols,1) = qfx(i,j) ! Surface constituent flux (kg/m2/s)
!Following variables are initialized to zero, they are the output from the "compute_eddy_diff" call
ustar8( :pcols) = 0.0_r8 ! Surface friction velocity (m/s)
pblh( :pcols) = 0.0_r8 ! Planetary boundary layer height (m )
ipbl( :pcols) = 0.0_r8 ! If 1, PBL is CL, while if 0, PBL is STL.
kpblh( :pcols) = 0.0_r8 ! Layer index containing PBL top within or at the base interface
wstarPBL(:pcols) = 0.0_r8 ! Convective velocity within PBL (m/s)
!!----------------------------------------------------------------------- !
!! Computation of eddy diffusivities - Select appropriate PBL scheme !
!!----------------------------------------------------------------------- !
select case (eddy_scheme)
case ( 'diag_TKE' )
!! ---------------------------------------------------------------- !
!! At first time step, have eddy_diff.F90:caleddy() use kvh=kvm=kvf !
!! This has to be done in compute_eddy_diff after kvf is calculated !
!! ---------------------------------------------------------------- !
kvinit = .false.
if(is_first_step) then
kvinit = .true.
endif
!! ---------------------------------------------- !
!! Get LW radiative heating out of physics buffer !
!! ---------------------------------------------- !
!! Retrieve eddy diffusivities for heat and momentum from physics buffer
!! from last timestep ( if first timestep, has been initialized by inidat.F90 )
kvm_in(:ncol,:) = kvm(:ncol,:)
kvh_in(:ncol,:) = kvh(:ncol,:)
call compute_eddy_diff( lchnk , pcols , pver , ncol , t8 , cloud(:,:,1) , ztodt , &
cloud(:,:,2), cloud(:,:,3), s8 , rpdel8 , cldn8 , qrl8 , wsedl8 , zm8 , zi8 , &
pmid8 , pint8 , u8 , v8 , taux , tauy , shflx , cflx(:,1), wstarent, &
nturb , ustar8 , pblh , kvm_in , kvh_in , kvm , kvh , kvq , cgh , &
cgs , tpert , qpert , wpert , tke8 , bprod , sprod , sfi , fqsatd , &
kvinit , tauresx , tauresy, ksrftms, ipbl(:), kpblh(:), wstarPBL(:) , turbtype , smaw )
!! ----------------------------------------------- !
!! Store TKE in pbuf for use by shallow convection !
!! ----------------------------------------------- !
!! Store updated kvh, kvm in pbuf to use here on the next timestep
do k = kts,kte+1
kflip = kte - k + 2
kvh3d(i,k,j) = kvh(1,kflip)
kvm3d(i,k,j) = kvm(1,kflip)
end do
end select
!!------------------------------------ !
!! Application of diffusivities !
!!------------------------------------ !
cloudtnd( :ncol,:,:) = cloud(:ncol,:,:)
stnd( :ncol,: ) = s8( :ncol,: )
wind_tends(:ncol,:,1) = u8( :ncol,: )
wind_tends(:ncol,:,2) = v8( :ncol,: )
!!------------------------------------------------------ !
!! Write profile output before applying diffusion scheme !
!!------------------------------------------------------ !
sl_prePBL(:ncol,:pver) = stnd(:ncol,:pver) - latvap * cloudtnd(:ncol,:pver,ixcldliq) &
- ( latvap + latice) * cloudtnd(:ncol,:pver,ixcldice)
qt_prePBL(:ncol,:pver) = cloudtnd(:ncol,:pver,1) + cloudtnd(:ncol,:pver,ixcldliq) &
+ cloudtnd(:ncol,:pver,ixcldice)
call aqsat( t8, pmid8, tem2, ftem, pcols, ncol, pver, 1, pver )
ftem_prePBL(:ncol,:) = cloud(:ncol,:,1)/ftem(:ncol,:)*100._r8
!! --------------------------------------------------------------------------------- !
!! Call the diffusivity solver and solve diffusion equation !
!! The final two arguments are optional function references to !
!! constituent-independent and constituent-dependent moleculuar diffusivity routines !
!! --------------------------------------------------------------------------------- !
!! Modification : We may need to output 'tautotx_im,tautoty_im' from below 'compute_vdiff' and
!! separately print out as diagnostic output, because these are different from
!! the explicit 'tautotx, tautoty' computed above.
!! Note that the output 'tauresx,tauresy' from below subroutines are fully implicit ones.
if( any(fieldlist_wet) ) then
call compute_vdiff( lchnk, pcols, pver, pcnst, ncol, pmid8, pint8, rpdel8, t8, ztodt, &
taux, tauy, shflx, cflx, ntop, nbot, kvh, kvm, kvq, cgs, cgh, zi8, ksrftms, qmincg, &
fieldlist_wet, wind_tends(:,:,1), wind_tends(:,:,2), cloudtnd, stnd, tautmsx, &
tautmsy, dtk, topflx, errstring, tauresx, tauresy, 1, do_molec_diff, &
compute_molec_diff, vd_lu_qdecomp)
end if
if( errstring .ne. '' ) then
call wrf_error_fatal(errstring)
endif
if( any( fieldlist_dry ) ) then
if( do_molec_diff ) then
errstring = "Design flaw: dry vdiff not currently supported with molecular diffusion"
call wrf_error_fatal(errstring)
end if
call compute_vdiff( lchnk, pcols, pver, pcnst, ncol, pmiddry8, pintdry8, rpdeldry8, t8, &
ztodt, taux, tauy, shflx, cflx, ntop, nbot, kvh, kvm, kvq, cgs, cgh, zi8, ksrftms, &
qmincg, fieldlist_dry, wind_tends(:,:,1), wind_tends(:,:,2), cloudtnd, stnd, tautmsx, &
tautmsy, dtk, topflx, errstring, tauresx, tauresy, 1, do_molec_diff, &
compute_molec_diff, vd_lu_qdecomp)
if( errstring .ne. '' ) call wrf_error_fatal(errstring)
end if
!! Store updated tauresx, tauresy to use here on the next timestep
tauresx2d(i,j) = tauresx(ncol)
tauresy2d(i,j) = tauresy(ncol)
#ifdef MODAL_AERO
!! Add the explicit surface fluxes to the lowest layer
!! Modification : I should check whether this explicit adding is consistent with
!! the treatment of other tracers.
!The following code is commented out as the diffusion for the Aerosols and other species is handled by CAMMGMP and WRF_CHEM's dry_dep_driver subroutines
!tmp1(:ncol) = ztodt * gravit * rpdel8(:ncol,pver)
!do m = 1, ntot_amode
! l = numptr_amode(m)
! cloudtnd(:ncol,pver,l) = cloudtnd(:ncol,pver,l) + tmp1(:ncol) * cflx(:ncol,l)
! do lspec = 1, nspec_amode(m)
! l = lmassptr_amode(lspec,m)
! cloudtnd(:ncol,pver,l) = cloudtnd(:ncol,pver,l) + tmp1(:ncol) * cflx(:ncol,l)
! enddo
!enddo
#endif
!! -------------------------------------------------------- !
!! Diagnostics and output writing after applying PBL scheme !
!! -------------------------------------------------------- !
sl(:ncol,:pver) = stnd(:ncol,:pver) - latvap * cloudtnd(:ncol,:pver,ixcldliq) &
- ( latvap + latice) * cloudtnd(:ncol,:pver,ixcldice)
qt(:ncol,:pver) = cloudtnd(:ncol,:pver,1) + cloudtnd(:ncol,:pver,ixcldliq) &
+ cloudtnd(:ncol,:pver,ixcldice)
!! --------------------------------------------------------------- !
!! Convert the new profiles into vertical diffusion tendencies. !
!! Convert KE dissipative heat change into "temperature" tendency. !
!! --------------------------------------------------------------- !
slten(:ncol,:) = ( sl(:ncol,:) - sl_prePBL(:ncol,:) ) * rztodt
qtten(:ncol,:) = ( qt(:ncol,:) - qt_prePBL(:ncol,:) ) * rztodt
stnd(:ncol,:) = ( stnd(:ncol,:) - s8(:ncol,:) ) * rztodt
wind_tends(:ncol,:,1) = ( wind_tends(:ncol,:,1) - u8(:ncol,:) ) * rztodt
wind_tends(:ncol,:,2) = ( wind_tends(:ncol,:,2) - v8(:ncol,:) ) * rztodt
cloudtnd(:ncol,:pver,:) = ( cloudtnd(:ncol,:pver,:) - cloud(:ncol,:pver,:) ) * rztodt
!! ----------------------------------------------------------- !
!! In order to perform 'pseudo-conservative varible diffusion' !
!! perform the following two stages: !
!! !
!! I. Re-set (1) 'qvten' by 'qtten', and 'qlten = qiten = 0' !
!! (2) 'sten' by 'slten', and !
!! (3) 'qlten = qiten = 0' !
!! !
!! II. Apply 'positive_moisture' !
!! !
!! ----------------------------------------------------------- !
if( eddy_scheme .eq. 'diag_TKE' .and. do_pseudocon_diff ) then
cloudtnd(:ncol,:pver,1) = qtten(:ncol,:pver)
stnd(:ncol,:pver) = slten(:ncol,:pver)
cloudtnd(:ncol,:pver,ixcldliq) = 0._r8
cloudtnd(:ncol,:pver,ixcldice) = 0._r8
cloudtnd(:ncol,:pver,ixnumliq) = 0._r8
cloudtnd(:ncol,:pver,ixnumice) = 0._r8
do ips = 1, ncol
do k = 1, pver
qv_pro(ips,k) = cloud(ips,k,1) + cloudtnd(ips,k,1) * ztodt
ql_pro(ips,k) = cloud(ips,k,ixcldliq) + cloudtnd(ips,k,ixcldliq) * ztodt
qi_pro(ips,k) = cloud(ips,k,ixcldice) + cloudtnd(ips,k,ixcldice) * ztodt
s_pro(ips,k) = s8(ips,k) + stnd(ips,k) * ztodt
t_pro(ips,k) = t8(ips,k) + (1._r8/cpair)*stnd(ips,k) * ztodt
end do
end do
call positive_moisture( cpair, latvap, latvap+latice, ncol, pver, ztodt, qmin(1), qmin(2), qmin(3), &
pdel8(:ncol,pver:1:-1), qv_pro(:ncol,pver:1:-1), ql_pro(:ncol,pver:1:-1), &
qi_pro(:ncol,pver:1:-1), t_pro(:ncol,pver:1:-1), s_pro(:ncol,pver:1:-1), &
cloudtnd(:ncol,pver:1:-1,1), cloudtnd(:ncol,pver:1:-1,ixcldliq), &
cloudtnd(:ncol,pver:1:-1,ixcldice), stnd(:ncol,pver:1:-1) )
end if
!! ----------------------------------------------------------------- !
!! Re-calculate diagnostic output variables after vertical diffusion !
!! ----------------------------------------------------------------- !
qv_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,1) + cloudtnd(:ncol,:pver,1) * ztodt
ql_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,ixcldliq) + cloudtnd(:ncol,:pver,ixcldliq) * ztodt
qi_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,ixcldice) + cloudtnd(:ncol,:pver,ixcldice) * ztodt
s_aft_PBL(:ncol,:pver) = s8(:ncol,:pver) + stnd(:ncol,:pver) * ztodt
t_aftPBL(:ncol,:pver) = ( s_aft_PBL(:ncol,:pver) - gravit*zm8(:ncol,:pver) ) / cpair
u_aft_PBL(:ncol,:pver) = u8(:ncol,:pver) + wind_tends(:ncol,:pver,1) * ztodt
v_aft_PBL(:ncol,:pver) = v8(:ncol,:pver) + wind_tends(:ncol,:pver,2) * ztodt
call aqsat( t_aftPBL, pmid8, tem2, ftem, pcols, ncol, pver, 1, pver )
ftem_aftPBL(:ncol,:pver) = qv_aft_PBL(:ncol,:pver) / ftem(:ncol,:pver) * 100._r8
tten(:ncol,:pver) = ( t_aftPBL(:ncol,:pver) - t8(:ncol,:pver) ) * rztodt
rhten(:ncol,:pver) = ( ftem_aftPBL(:ncol,:pver) - ftem_prePBL(:ncol,:pver) ) * rztodt
!Post processing of the output from CAM's parameterization
do k=kts,kte
kflip = kte-k+1
rublten(i,k,j) = wind_tends(1,kflip,1)
rvblten(i,k,j) = wind_tends(1,kflip,2)
rthblten(i,k,j) = stnd(1,kflip)/cpair/exner8(1,kflip)
multFrc = 1._r8 + qv_curr(i,k,j)
rqvblten(i,k,j) = cloudtnd(1,kflip,1 ) * multFrc * multFrc
rqcblten(i,k,j) = cloudtnd(1,kflip,ixcldliq) * multFrc
rqiblten(i,k,j) = cloudtnd(1,kflip,ixcldice) * multFrc
!*Important* : ixnumliq is mixed in the dropmixnuc, therefore ixnumliq is NOT mixed here (ONLY if CAMMGMP is used for mp_physics)
rqniblten(i,k,j) = cloudtnd(1,kflip,ixnumice) * multFrc
!Diffusivity coefficients at the midpoints
kp1 = k + 1
end do
do k = kts,kte+1
kflip = kte - k + 2
tke_pbl(i,k,j) = tke8(1,kflip) !TKE at the interfaces
turbtype3d(i,k,j) = turbtype(1,kflip)
smaw3d(i,k,j) = smaw(1,kflip)
end do
kpbl2d(i,j) = kte - int(kpblh(1)) + 1 !int(kpblh(1))
pblh2d(i,j) = pblh( 1)
tpert2d(i,j) = tpert(1)
qpert2d(i,j) = qpert(1)
wpert2d(i,j) = wpert(1)
!End Post processing of the output from CAM
enddo !loop of i
enddo !loop of j
return
end subroutine camuwpbl
!-----------------------------------------------------------------------------------------
subroutine positive_moisture( cp, xlv, xls, ncol, mkx, dt, qvmin, qlmin, qimin, &
dp, qv, ql, qi, t, s, qvten, qlten, qiten, sten )
!! ------------------------------------------------------------------------------- !
!! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, !
!! force them to be larger than minimum value by (1) condensating water vapor !
!! into liquid or ice, and (2) by transporting water vapor from the very lower !
!! layer. '2._r8' is multiplied to the minimum values for safety. !
!! Update final state variables and tendencies associated with this correction. !
!! If any condensation happens, update (s,t) too. !
!! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
!! input tendencies. !
!! Be careful the order of k : '1': near-surface layer, 'mkx' : top layer !
!! ------------------------------------------------------------------------------- !
!-----------------------------------------------------------------------------------------
implicit none
integer, intent(in) :: ncol, mkx
real(r8), intent(in) :: cp, xlv, xls
real(r8), intent(in) :: dt, qvmin, qlmin, qimin
real(r8), intent(in) :: dp(ncol,mkx)
real(r8), intent(inout) :: qv(ncol,mkx), ql(ncol,mkx), qi(ncol,mkx), t(ncol,mkx), s(ncol,mkx)
real(r8), intent(inout) :: qvten(ncol,mkx), qlten(ncol,mkx), qiten(ncol,mkx), sten(ncol,mkx)
integer i, k
real(r8) dql, dqi, dqv, sum, aa, dum
!! Modification : I should check whether this is exactly same as the one used in
!! shallow convection and cloud macrophysics.
do i = 1, ncol
do k = mkx, 1, -1 ! From the top to the 1st (lowest) layer from the surface
dql = max(0._r8,1._r8*qlmin-ql(i,k))
dqi = max(0._r8,1._r8*qimin-qi(i,k))
qlten(i,k) = qlten(i,k) + dql/dt
qiten(i,k) = qiten(i,k) + dqi/dt
qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
sten(i,k) = sten(i,k) + xlv * (dql/dt) + xls * (dqi/dt)
ql(i,k) = ql(i,k) + dql
qi(i,k) = qi(i,k) + dqi
qv(i,k) = qv(i,k) - dql - dqi
s(i,k) = s(i,k) + xlv * dql + xls * dqi
t(i,k) = t(i,k) + (xlv * dql + xls * dqi)/cp
dqv = max(0._r8,1._r8*qvmin-qv(i,k))
qvten(i,k) = qvten(i,k) + dqv/dt
qv(i,k) = qv(i,k) + dqv
if( k .ne. 1 ) then
qv(i,k-1) = qv(i,k-1) - dqv*dp(i,k)/dp(i,k-1)
qvten(i,k-1) = qvten(i,k-1) - dqv*dp(i,k)/dp(i,k-1)/dt
endif
qv(i,k) = max(qv(i,k),qvmin)
ql(i,k) = max(ql(i,k),qlmin)
qi(i,k) = max(qi(i,k),qimin)
end do
!! Extra moisture used to satisfy 'qv(i,1)=qvmin' is proportionally
!! extracted from all the layers that has 'qv > 2*qvmin'. This fully
!! preserves column moisture.
if( dqv .gt. 1.e-20_r8 ) then
sum = 0._r8
do k = 1, mkx
if( qv(i,k) .gt. 2._r8*qvmin ) sum = sum + qv(i,k)*dp(i,k)
enddo
aa = dqv*dp(i,1)/max(1.e-20_r8,sum)
if( aa .lt. 0.5_r8 ) then
do k = 1, mkx
if( qv(i,k) .gt. 2._r8*qvmin ) then
dum = aa*qv(i,k)
qv(i,k) = qv(i,k) - dum
qvten(i,k) = qvten(i,k) - dum/dt
endif
enddo
else
write(iulog,*) 'Full positive_moisture is impossible in vertical_diffusion'
call wrf_message(iulog)
endif
endif
end do
return
end subroutine positive_moisture
!-----------------------------------------------------------------------------------------
subroutine camuwpblinit(rublten,rvblten,rthblten,rqvblten, &
restart,tke_pbl,is_CAMMGMP_used, &
ids,ide,jds,jde,kds,kde, &
ims,ime,jms,jme,kms,kme, &
its,ite,jts,jte,kts,kte)
!!------------------------------------------------------------------!
!! Initialization of time independent fields for vertical diffusion !
!! Calls initialization routines for subsidiary modules !
!!----------------------------------------------------------------- !
!This subroutine is based on vertical_diffusion_init of CAM. This subroutine
!initializes variables for vertical diffusion subroutine calls. The layout
!is kept similar to the original CAM subroutine to facillitate future adaptations.
!Called by : physics_init.F
!Calls: vd_register, cnst_get_ind, wrf_error_fatal,init_eddy_diff,init_tms,init_vdiff
!-----------------------------------------------------------------------------------------
use eddy_diff, only : init_eddy_diff
use molec_diff, only : init_molec_diff
use trb_mtn_stress, only : init_tms
use diffusion_solver, only : init_vdiff, vdiff_select
use constituents, only : cnst_get_ind, cnst_get_type_byind, cnst_name
use module_cam_support, only : masterproc
use module_model_constants, only : epsq2
#ifdef MODAL_AERO
use modal_aero_data
#endif
implicit none
!-------------------------------------------------------------------------------------!
!Input and output variables !
!-------------------------------------------------------------------------------------!
logical,intent(in) :: restart,is_CAMMGMP_used
integer,intent(in) :: ids,ide,jds,jde,kds,kde
integer,intent(in) :: ims,ime,jms,jme,kms,kme
integer,intent(in) :: its,ite,jts,jte,kts,kte
real,dimension(ims:ime,kms:kme,jms:jme),intent(inout) :: &
rublten, rvblten, rthblten, rqvblten
real,dimension(ims:ime,kms:kme,jms:jme),intent(out) :: TKE_PBL
!-------------------------------------------------------------------------------------!
!Local Variables !
!-------------------------------------------------------------------------------------!
integer :: i,j,k,itf,jtf,ktf
integer :: ntop_eddy !! Top interface level to which eddy vertical diffusion is applied ( = 1 )
integer :: nbot_eddy !! Bottom interface level to which eddy vertical diffusion is applied ( = pver )
integer :: ntop_molec !! Top interface level to which molecular vertical diffusion is applied ( = 1 )
integer :: nbot_molec !! Bottom interface level to which molecular vertical diffusion is applied
character(128) :: errstring !! Error status for init_vdiff
real(r8) :: hypm(kte) !! reference state midpoint pressures
#ifdef MODAL_AERO
integer :: m, l
#endif
jtf = min(jte,jde-1)
ktf = min(kte,kde-1)
itf = min(ite,ide-1)
!Map CAM veritcal level variables
pver = ktf - kts + 1
pverp = pver + 1
!Initialize flags and add constituents
call vd_register()
!! ----------------------------------------------------------------- !
!! Get indices of cloud liquid and ice within the constituents array !
!! ----------------------------------------------------------------- !
call cnst_get_ind( 'CLDLIQ', ixcldliq )
call cnst_get_ind( 'CLDICE', ixcldice )
if( microp_scheme .eq. 'MG' ) then
call cnst_get_ind( 'NUMLIQ', ixnumliq )
call cnst_get_ind( 'NUMICE', ixnumice )
endif
if (masterproc) then
write(iulog,*)'Initializing vertical diffusion (vertical_diffusion_init)'
call wrf_debug(1,iulog)
end if
!! ---------------------------------------------------------------------------------------- !
!! Initialize molecular diffusivity module !
!! Molecular diffusion turned on above ~60 km (50 Pa) if model top is above ~90 km (.1 Pa). !
!! Note that computing molecular diffusivities is a trivial expense, but constituent !
!! diffusivities depend on their molecular weights. Decomposing the diffusion matric !
!! for each constituent is a needless expense unless the diffusivity is significant. !
!! ---------------------------------------------------------------------------------------- !
ntop_molec = 1 !! Should always be 1
nbot_molec = 0 !! Should be set below about 70 km
!! ---------------------------------- !
!! Initialize eddy diffusivity module !
!! ---------------------------------- !
ntop_eddy = 1 !! No reason not to make this 1, if > 1, must be <= nbot_molec
nbot_eddy = pver !! Should always be pver
if( masterproc ) write(iulog,fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY =', ntop_eddy, 'NBOT_EDDY =', nbot_eddy
call wrf_debug(1,iulog)
select case ( eddy_scheme )
case ( 'diag_TKE' )
if( masterproc ) write(iulog,*) 'vertical_diffusion_init: eddy_diffusivity scheme'
call wrf_debug(1,iulog)
if( masterproc ) write(iulog,*) 'UW Moist Turbulence Scheme by Bretherton and Park'
call wrf_debug(1,iulog)
!! Check compatibility of eddy and shallow scheme
if( shallow_scheme .ne. 'UW' ) then
write(iulog,*) 'ERROR: shallow convection scheme ', shallow_scheme,' is incompatible with eddy scheme ', eddy_scheme
call wrf_message(iulog)
call wrf_error_fatal( 'convect_shallow_init: shallow_scheme and eddy_scheme are incompatible' )
endif
call init_eddy_diff( r8, pver, gravit, cpair, rair, zvir, latvap, latice, &
ntop_eddy, nbot_eddy, hypm, karman )
if( masterproc ) write(iulog,*) 'vertical_diffusion: nturb, ntop_eddy, nbot_eddy ', nturb, ntop_eddy, nbot_eddy
call wrf_debug(1,iulog)
end select
!!-------------------------------------------------------------------------------------!
!! The vertical diffusion solver must operate !
!! over the full range of molecular and eddy diffusion !
!!-------------------------------------------------------------------------------------!
ntop = min(ntop_molec,ntop_eddy)
nbot = max(nbot_molec,nbot_eddy)
!! ------------------------------------------- !
!! Initialize turbulent mountain stress module !
!! ------------------------------------------- !
if( do_tms ) then
call init_tms( r8, tms_orocnst, tms_z0fac, karman, gravit, rair )
if (masterproc) then
write(iulog,*)'Using turbulent mountain stress module'
call wrf_message(iulog)
write(iulog,*)' tms_orocnst = ',tms_orocnst
call wrf_message(iulog)
end if
endif
!! ---------------------------------- !
!! Initialize diffusion solver module !
!! ---------------------------------- !
call init_vdiff( r8, pcnst, rair, gravit, fieldlist_wet, fieldlist_dry, errstring )