wrf-fire /wrfv2_fire/phys/module_bl_acm.F

Language Fortran 77 Lines 1211
MD5 Hash 596dc859b7d9c664a4454fbd0a5733c4 Estimated Cost $19,554 (why?)
Repository git://github.com/jbeezley/wrf-fire.git View Raw File View Project SPDX
   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
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_bl_acm

!  USE module_model_constants

  REAL, PARAMETER      :: RIC    = 0.25                ! critical Richardson number
  REAL, PARAMETER      :: CRANKP = 0.5                 ! CRANK-NIC PARAMETER

CONTAINS

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   SUBROUTINE ACMPBL(XTIME,    DTPBL,    ZNW,   SIGMAH,               &
                     U3D,      V3D,      PP3D,  DZ8W, TH3D, T3D,      &
                     QV3D,     QC3D,     QI3D,  RR3D,                 &
                     UST,      HFX,      QFX,   TSK,                  &
                     PSFC,     EP1,      G,                           &
                     ROVCP,    RD,       CPD,                         &
                     PBLH,     KPBL2D,   REGIME,                      &
                     GZ1OZ0,   WSPD,     PSIM, MUT,                   &
                     RUBLTEN,  RVBLTEN,  RTHBLTEN,                    &
                     RQVBLTEN, RQCBLTEN, RQIBLTEN,                    &
                     ids,ide, jds,jde, kds,kde,                       &
                     ims,ime, jms,jme, kms,kme,                       &
                     its,ite, jts,jte, kts,kte)
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

!   THIS MODULE COMPUTES VERTICAL MIXING IN AND ABOVE THE PBL ACCORDING TO 
!   THE ASYMMETRICAL CONVECTIVE MODEL, VERSION 2  (ACM2), WHICH IS A COMBINED 
!   LOCAL NON-LOCAL CLOSURE SCHEME BASED ON THE ORIGINAL ACM (PLEIM AND CHANG 1992)
!
!   REFERENCES: 
!   Pleim (2007) A combined local and non-local closure model for the atmospheric
!                boundary layer.  Part1: Model description and testing.  
!                JAMC, 46, 1383-1395
!   Pleim (2007) A combined local and non-local closure model for the atmospheric
!                boundary layer.  Part2: Application and evaluation in a mesoscale
!                meteorology model.  JAMC, 46, 1396-1409
!
!  REVISION HISTORY:
!     AX        3/2005   - developed WRF version based on the MM5 PX LSM
!     RG and JP 7/2006   - Finished WRF adaptation
!     JP 12/2011 12/2011 - ACM2 modified so it's not dependent on first layer thickness. 
!
!**********************************************************************
!   ARGUMENT LIST:
!
!... Inputs:
!-- XTIME           Time since simulation start (min)
!-- DTPBL           PBL time step
!-- ZNW             Sigma at full layer
!-- SIGMAH          Sigma at half layer
!-- U3D             3D u-velocity interpolated to theta points (m/s)
!-- V3D             3D v-velocity interpolated to theta points (m/s)
!-- PP3D            Pressure at half levels (Pa)
!-- DZ8W            dz between full levels (m)
!-- TH3D            Potential Temperature (K)
!-- T3D             Temperature (K)
!-- QV3D            3D water vapor mixing ratio (Kg/Kg)
!-- QC3D            3D cloud mixing ratio (Kg/Kg)
!-- QI3D            3D ice mixing ratio (Kg/Kg)
!-- RR3D            3D dry air density (kg/m^3)
!-- UST             Friction Velocity (m/s)
!-- HFX		    Upward heat flux at the surface (w/m^2)
!-- QFX		    Upward moisture flux at the surface (Kg/m^2/s)
!-- TSK             Surface temperature (K)
!-- PSFC            Pressure at the surface (Pa)
!-- EP1             Constant for virtual temperature (r_v/r_d-1) (dimensionless)
!-- G               Gravity (m/s^2)
!-- ROVCP           r/cp
!-- RD              gas constant for dry air (j/kg/k)
!-- CPD             heat capacity at constant pressure for dry air (j/kg/k)
!-- GZ1OZ0          log(z/z0) where z0 is roughness length
!-- WSPD            wind speed at lowest model level (m/s)
!-- PSIM            similarity stability function for momentum
!-- MUT             Total Mu : Psfc - Ptop
 
!-- ids             start index for i in domain
!-- ide             end index for i in domain
!-- jds             start index for j in domain
!-- jde             end index for j in domain
!-- kds             start index for k in domain
!-- kde             end index for k in domain
!-- ims             start index for i in memory
!-- ime             end index for i in memory
!-- jms             start index for j in memory
!-- jme             end index for j in memory
!-- kms             start index for k in memory
!-- kme             end index for k in memory
!-- jts             start index for j in tile
!-- jte             end index for j in tile
!-- kts             start index for k in tile
!-- kte             end index for k in tile
!
!... Outputs: 
!-- PBLH            PBL height (m)
!-- KPBL2D          K index for PBL layer
!-- REGIME          Flag indicating PBL regime (stable, unstable, etc.)
!-- RUBLTEN         U tendency due to PBL parameterization (m/s^2)
!-- RVBLTEN         V tendency due to PBL parameterization (m/s^2)
!-- RTHBLTEN        Theta tendency due to PBL parameterization (K/s)
!-- RQVBLTEN        Qv tendency due to PBL parameterization (kg/kg/s)
!-- RQCBLTEN        Qc tendency due to PBL parameterization (kg/kg/s)
!-- RQIBLTEN        Qi tendency due to PBL parameterization (kg/kg/s)
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
     IMPLICIT NONE

!.......Arguments
! DECLARATIONS - INTEGER
    INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
                                      ims,ime, jms,jme, kms,kme, &
                                      its,ite, jts,jte, kts,kte, XTIME

! DECLARATIONS - REAL
    REAL,                                INTENT(IN)  ::  DTPBL, EP1,   &
                                                        G, ROVCP, RD, CPD

    REAL,    DIMENSION( kms:kme ),       INTENT(IN)  :: ZNW, SIGMAH

    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ),                         &
             INTENT(IN) ::                              U3D, V3D,            &
                                                        PP3D, DZ8W, T3D,     &
                                                        QV3D, QC3D, QI3D,    &
                                                        RR3D, TH3D

    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: PSIM, GZ1OZ0,     &
                                                          HFX, QFX, TSK,    &
                                                          PSFC, WSPD, MUT

    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::  PBLH, REGIME, UST

    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ),                         &
             INTENT(INOUT)   ::                         RUBLTEN, RVBLTEN,    &
                                                        RTHBLTEN, RQVBLTEN,  &
                                                        RQCBLTEN, RQIBLTEN

    INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(OUT  ) ::  KPBL2D

!... Local Variables

!... Integer
      INTEGER :: J, K
!... Real
      REAL, DIMENSION( kts:kte ) :: DSIGH, DSIGHI, DSIGFI
      REAL, DIMENSION( 0:kte )   :: SIGMAF
      REAL  RDT
      REAL, PARAMETER :: KARMAN = 0.4
!...

   RDT = 1.0 / DTPBL

   DO K = 1, kte
     SIGMAF(K-1) = ZNW(K)
   ENDDO
   SIGMAF(kte) = 0.0

   DO K = 1, kte
     DSIGH(K)  = SIGMAF(K) - SIGMAF(K-1)
     DSIGHI(K) = 1.0 / DSIGH(K)
   ENDDO

   DO K = kts,kte-1
     DSIGFI(K) = 1.0 / (SIGMAH(K+1) - SIGMAH(K))
   ENDDO

   DSIGFI(kte) = DSIGFI(kte-1)
   
   DO j = jts,jte
      CALL ACM2D(j=J,xtime=XTIME, dtpbl=DTPBL, sigmaf=SIGMAF, sigmah=SIGMAH    &
              ,dsigfi=DSIGFI,dsighi=DSIGHI,dsigh=DSIGH             &
              ,us=u3d(ims,kms,j),vs=v3d(ims,kms,j)                 &
              ,theta=th3d(ims,kms,j),tt=t3d(ims,kms,j)             &
              ,qvs=qv3d(ims,kms,j),qcs=qc3d(ims,kms,j)             &
              ,qis=qi3d(ims,kms,j),dzf=DZ8W(ims,kms,j)             &
              ,densx=RR3D(ims,kms,j)                               &
              ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j)     &
              ,ttnp=rthblten(ims,kms,j),qvtnp=rqvblten(ims,kms,j)  &
              ,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &
              ,cpd=cpd,g=g,rovcp=rovcp,rd=rd,rdt=rdt               &
              ,psfcpa=psfc(ims,j),ust=ust(ims,j)                   &
              ,pbl=pblh(ims,j)                                     &
              ,regime=regime(ims,j),psim=psim(ims,j)               &
              ,hfx=hfx(ims,j),qfx=qfx(ims,j)                       &
              ,tg=tsk(ims,j),gz1oz0=gz1oz0(ims,j)                  &
              ,wspd=wspd(ims,j) ,klpbl=kpbl2d(ims,j)               &
              ,mut=mut(ims,j)                                      &
              ,ep1=ep1,karman=karman                               &
              ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde   &
              ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme   &
              ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
   ENDDO

   END SUBROUTINE ACMPBL
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------


!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   SUBROUTINE ACM2D(j,XTIME, DTPBL, sigmaf, sigmah          &
              ,dsigfi,dsighi,dsigh                          &
              ,us,vs,theta,tt,qvs,qcs,qis                   &
              ,dzf,densx,utnp,vtnp,ttnp,qvtnp,qctnp,qitnp   &
              ,cpd,g,rovcp,rd,rdt,psfcpa,ust                &
              ,pbl,regime,psim                              &
              ,hfx,qfx,tg,gz1oz0,wspd ,klpbl                &
              ,mut                                          &
              ,ep1,karman                                   &
              ,ids,ide, jds,jde, kds,kde   &
              ,ims,ime, jms,jme, kms,kme   &
              ,its,ite, jts,jte, kts,kte   )

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
      IMPLICIT NONE

!.......Arguments

!... Real
      REAL, DIMENSION( 0:kte ),             INTENT(IN)  :: SIGMAF
      REAL, DIMENSION( kms:kme ),           INTENT(IN)  :: SIGMAH
      REAL, DIMENSION( kts:kte ),           INTENT(IN)  :: DSIGH, DSIGHI, DSIGFI
      REAL ,                                INTENT(IN)  :: DTPBL, G, RD,ep1,karman,CPD,ROVCP,RDT
      REAL , DIMENSION( ims:ime ),          INTENT(INOUT)  :: PBL, UST
      
      REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, THETA, TT,   &
                                                           QVS, QCS, QIS, DENSX
      REAL,  DIMENSION( ims:ime, kms:kme ), intent(in)  :: DZF
      REAL,  DIMENSION( ims:ime, kms:kme ), intent(inout)   ::  utnp, &
							        vtnp, &
							        ttnp, &
							        qvtnp, &
							        qctnp, &
							        qitnp
      real,     dimension( ims:ime ), intent(in   )   ::   psfcpa
      real,     dimension( ims:ime ), intent(in   )   ::   tg
      real,     dimension( ims:ime ), intent(inout)   ::   regime
      real,     dimension( ims:ime ), intent(in)      ::   wspd, psim, gz1oz0
      real,     dimension( ims:ime ), intent(in)      ::   hfx, qfx
      real,     dimension( ims:ime ), intent(in)      ::   mut

!... Integer
      INTEGER, DIMENSION( ims:ime ),       INTENT(OUT):: KLPBL
      INTEGER,  INTENT(IN)      ::      XTIME
      integer,  intent(in   )   ::      ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, j
!--------------------------------------------------------------------
!--Local 
      INTEGER I, K     
      INTEGER :: KPBLHT
      INTEGER, DIMENSION( its:ite ) :: KPBLH, NOCONV

!... Real
      REAL    ::  TVCON, WSS, TCONV, TH1, TOG, DTMP, WSSQ
      REAL    ::  psix, THV1
      REAL, DIMENSION( its:ite )          :: FINT, PSTAR, CPAIR
      REAL, DIMENSION( its:ite, kts:kte ) :: THETAV, RIB,             &
                                             EDDYZ, UX, VX, THETAX,   &
                                             QVX, QCX, QIX, ZA
      REAL, DIMENSION( its:ite, 0:kte )   :: ZF
      REAL,    DIMENSION( its:ite)           :: WST, TST, QST, USTM, TSTV
      REAL,    DIMENSION( its:ite )          :: PBLSIG, MOL
      REAL    ::  FINTT, ZMIX, UMIX, VMIX
      REAL    ::  TMPFX, TMPVTCON, TMPP, TMPTHS, TMPTH1, TMPVCONV, WS1, DTH
	REAL    ::  A,TST12,RL,ZFUNC
!    REAL, PARAMETER :: KARMAN = 0.4

!... Integer
      INTEGER :: KL, jtf, ktf, itf, KMIX, KSRC
!...
        character*256 :: message
!-----initialize vertical tendencies and

      DO i = its,ite
        DO k = kts,kte
          utnp(i,k) = 0.
          vtnp(i,k) = 0.
          ttnp(i,k) = 0.
        ENDDO
      ENDDO

      DO k = kts,kte
        DO i = its,ite
          qvtnp(i,k) = 0.
        ENDDO
      ENDDO

      DO k = kts,kte
        DO i = its,ite
          qctnp(i,k) = 0.
          qitnp(i,k) = 0.
        ENDDO
      ENDDO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!  Compute Micromet Scaling variables, not availiable in WRF for ACM
     DO I = its,ite
           CPAIR(I)  = CPD * (1.0 + 0.84 * QVS(I,1))                    ! J/(K KG)
           TMPFX     = HFX(I)  / (cpair(i) * DENSX(I,1))
           TMPVTCON  = 1.0 + EP1 * QVS(I,1)                             ! COnversion factor for virtual temperature
           WS1       = SQRT(US(I,1)**2 + VS(I,1)**2)                    ! Level 1 wind speed
           TST(I)    = -TMPFX / UST(I)
           QST(I)    = -QFX(I) / (UST(I)*DENSX(I,1))
           USTM(I)   = UST(I) * WS1 / wspd(i)
           THV1      = TMPVTCON * THETA(I,1) 
           TSTV(I)   = TST(I)*TMPVTCON + THV1*EP1*QST(I)
           IF(ABS(TSTV(I)).LT.1.0E-6) THEN
             TSTV(I) = SIGN(1.0E-6,TSTV(I))
           ENDIF
           MOL(I)    = THV1 * UST(i)**2/(KARMAN*G*TSTV(I))
           WST(I)    = UST(I) * (PBL(I)/(KARMAN*ABS(MOL(I)))) ** 0.333333       
           PSTAR(I)  =  MUT(I)/1000.                                     ! P* in cb 
     ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!... Compute PBL height

!... compute the height of full- and half-sigma level above ground level
     DO I = its,ite
       ZF(I,0)    = 0.0
       KLPBL(I)   = 1
     ENDDO

     DO K = kts, kte
       DO I = its,ite
         ZF(I,K) = DZF(I,K) + ZF(I,K-1)
         ZA(I,K) = 0.5 * (ZF(I,K) + ZF(I,K-1))
       ENDDO
     ENDDO

     DO K = kts, kte
       DO I = its,ite
         TVCON       = 1.0 + EP1 * QVS(I,K)
         THETAV(I,K) = THETA(I,K) * TVCON
       ENDDO
     ENDDO


!...  COMPUTE PBL WHERE RICHARDSON NUMBER = RIC (0.25) HOLTSLAG ET AL 1990  
     DO 100 I = its,ite
       DO K = 1,kte
         KSRC = K
         IF (SIGMAF(K).lT.0.9955) GO TO 69
       ENDDO
69     CONTINUE
       TH1 = 0.0
       DO K = 1,KSRC
         TH1 = TH1 + THETAV(I,K)  
       ENDDO  
       TH1 = TH1/KSRC
       IF(MOL(I).LT.0.0 .AND. XTIME.GT.1) then
         WSS   = (UST(I) ** 3 + 0.6 * WST(I) ** 3) ** 0.33333
         TCONV = -8.5 * UST(I) * TSTV(I) / WSS
         TH1   = TH1 + TCONV
       ENDIF

99     KMIX = 1
       DO K = 1,kte
         DTMP   = THETAV(I,K) - TH1
         IF (DTMP.LT.0.0) KMIX = K
       ENDDO
       IF(KMIX.GT.1) THEN
         FINTT = (TH1 - THETAV(I,KMIX)) / (THETAV(I,KMIX+1)               &
               - THETAV(I,KMIX))
         ZMIX = FINTT * (ZA(I,KMIX+1)-ZA(I,KMIX)) + ZA(I,KMIX)
         UMIX = FINTT * (US(I,KMIX+1)-US(I,KMIX)) + US(I,KMIX)
         VMIX = FINTT * (VS(I,KMIX+1)-VS(I,KMIX)) + VS(I,KMIX)
       ELSE
         ZMIX = ZA(I,1)
         UMIX = US(I,1)
         VMIX = VS(I,1)
       ENDIF
       DO K = KMIX,kte
         DTMP   = THETAV(I,K) - TH1
         TOG = 0.5 * (THETAV(I,K) + TH1) / G
         WSSQ = (US(I,K)-UMIX)**2                                     &
              + (VS(I,K)-VMIX)**2
         IF (KMIX == 1) WSSQ = WSSQ + 100.*UST(I)*UST(I) 
         WSSQ = MAX( WSSQ, 0.1 )
         RIB(I,K) = ABS(ZA(I,K)-ZMIX) * DTMP / (TOG * WSSQ)
         IF (RIB(I,K) .GE. RIC) GO TO 201
       ENDDO

       write (message, *)' RIBX never exceeds RIC, RIB(i,kte) = ',rib(i,5),        &
               ' THETAV(i,1) = ',thetav(i,1),' MOL=',mol(i),            &
               ' TCONV = ',TCONV,' WST = ',WST(I),                      &
               ' KMIX = ',kmix,' UST = ',UST(I),                       &
               ' TST = ',TST(I),' U,V = ',US(I,1),VS(I,1),              &
               ' I,J=',I,J
       CALL wrf_error_fatal ( message )
201    CONTINUE

       KPBLH(I) = K

100  CONTINUE

     DO I = its,ite
       IF (KPBLH(I) .NE. 1) THEN
!---------INTERPOLATE BETWEEN LEVELS -- jp 7/93
         FINT(I) = (RIC - RIB(I,KPBLH(I)-1)) / (RIB(I,KPBLH(I)) -       &
                    RIB(I,KPBLH(I)-1))
         IF (FINT(I) .GT. 0.5) THEN
           KPBLHT  = KPBLH(I)
           FINT(I) = FINT(I) - 0.5
         ELSE
           KPBLHT  = KPBLH(I) - 1
           FINT(I) = FINT(I) + 0.5
         ENDIF
         PBL(I)  = FINT(I) * (ZF(I,KPBLHT) - ZF(I,KPBLHT-1)) +          &
                     ZF(I,KPBLHT-1)
         KLPBL(I) = KPBLHT
         PBLSIG(I)   = FINT(I) * DSIGH(KPBLHT) + SIGMAF(KPBLHT-1)    ! sigma at PBL height
       ELSE
         KLPBL(I) = 1
         PBL(I)    = ZF(I,1)                                                  
         PBLSIG(I)   = SIGMAF(1)                                             
       ENDIF

     ENDDO

     DO I = its,ite       
       NOCONV(I) = 0
       
! Check for CBL and identify conv. vs. non-conv cells
       IF (PBL(I) / MOL(I) .LT. -0.02 .AND. KLPBL(I) .GT. 3        &
           .AND. THETAV(I,1) .GT. THETAV(I,2) .AND. XTIME .GT. 1) THEN
          NOCONV(I)   = 1
          REGIME(I) = 4.0                     ! FREE CONVECTIVE - ACM
       ENDIF
     ENDDO

!... Calculate Kz
     CALL EDDYX(DTPBL, ZF,  ZA,     MOL, PBL,  UST,                &
                US,    VS,  TT,  THETAV, DENSX, PSTAR,              &
                QVS,   QCS, QIS, DSIGFI, G, RD, CPAIR,              &
                EDDYZ, its,ite, kts,kte,ims,ime, kms,kme)

     CALL ACM(DTPBL, PSTAR,  NOCONV, SIGMAF, DSIGH, DSIGHI, J,      &
                 KLPBL, PBL,   PBLSIG, MOL,  UST,                  &
                 TST, QST,  USTM,   EDDYZ,  DENSX,                  &
                 US,    VS,     THETA,  QVS,    QCS,    QIS,        &
                 UX,    VX,     THETAX, QVX,    QCX,    QIX,        &
                 ids,ide, jds,jde, kds,kde,                         &
                 ims,ime, jms,jme, kms,kme,                         &
                 its,ite, jts,jte, kts,kte)

!... Calculate tendency due to PBL parameterization

     DO K = kts, kte
       DO I = its, ite
         UTNP(I,K)  = UTNP(I,K) + (UX(I,K) - US(I,K)) * RDT
         VTNP(I,K)  = VTNP(I,K) + (VX(I,K) - VS(I,K)) * RDT
         TTNP(I,K)  = TTNP(I,K) + (THETAX(I,K) - THETA(I,K)) * RDT
         QVTNP(I,K) = QVTNP(I,K) + (QVX(I,K) - QVS(I,K)) * RDT
         QCTNP(I,K) = QCTNP(I,K) + (QCX(I,K) - QCS(I,K)) * RDT
         QITNP(I,K) = QITNP(I,K) + (QIX(I,K) - QIS(I,K)) * RDT
       ENDDO
     ENDDO

   END SUBROUTINE ACM2D
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   SUBROUTINE ACMINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,           &
                      RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR,       &
                      restart, allowed_to_read ,                   &
                      ids, ide, jds, jde, kds, kde,                &
                      ims, ime, jms, jme, kms, kme,                &
                      its, ite, jts, jte, kts, kte                 )
!-----------------------------------------------------------------------
!
!    This subroutine is for preparing ACM PBL variables. 
!    Called from module_physics_init.F
!
!  REVISION HISTORY:
!     AX     3/2005 - Originally developed
!-----------------------------------------------------------------------
!   ARGUMENT LIST:
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
     IMPLICIT NONE
!
   LOGICAL , INTENT(IN)          :: restart , allowed_to_read

   INTEGER , INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
                                     ims, ime, jms, jme, kms, kme, &
                                     its, ite, jts, jte, kts, kte

   INTEGER , INTENT(IN)          ::  P_QI,P_FIRST_SCALAR

!   REAL , DIMENSION( kms:kme ), INTENT(IN)  :: SHALF
   REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::         &
                                                         RUBLTEN, &
                                                         RVBLTEN, &
                                                         RTHBLTEN, &
                                                         RQVBLTEN, &
                                                         RQCBLTEN, & 
                                                         RQIBLTEN

!... Local Variables
   INTEGER :: i, j, k, itf, jtf, ktf

!
   jtf=min0(jte,jde-1)
   ktf=min0(kte,kde-1)
   itf=min0(ite,ide-1)

   IF(.not.restart)THEN
     DO j=jts,jtf
     DO k=kts,ktf
     DO i=its,itf
        RUBLTEN(i,k,j)=0.
        RVBLTEN(i,k,j)=0.
        RTHBLTEN(i,k,j)=0.
        RQVBLTEN(i,k,j)=0.
        RQCBLTEN(i,k,j)=0.
     ENDDO
     ENDDO
     ENDDO
   ENDIF

   IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
      DO j=jts,jtf
      DO k=kts,ktf
      DO i=its,itf
         RQIBLTEN(i,k,j)=0.
      ENDDO
      ENDDO
      ENDDO
   ENDIF


   END SUBROUTINE acminit
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!-------------------------------------------------------------------          
   SUBROUTINE EDDYX(DTPBL, ZF,  ZA,     MOL, PBL,  UST,               &
                    US,    VS,  TT,  THETAV, DENSX, PSTAR,            &
                    QVS,   QCS, QIS, DSIGFI, G, RD, CPAIR,            &
                    EDDYZ, its,ite, kts,kte,ims,ime,kms,kme )


!**********************************************************************
!   Two methods for computing Kz:
!   1.  Boundary scaling similar to Holtslag and Boville (1993)
!   2.  Local Kz computed as function of local Richardson # and vertical 
!       wind shear, similar to LIU & CARROLL (1996)
!
!**********************************************************************
!
!-- DTPBL           time step of the minor loop for the land-surface/pbl model
!-- ZF              height of full sigma level
!-- ZA              height of half sigma level
!-- MOL             Monin-Obukhov length in 1D form
!-- PBL             PBL height in 1D form
!-- UST             friction velocity U* in 1D form (m/s)
!-- US              U wind 
!-- VS              V wind
!-- TT              temperature
!-- THETAV          potential virtual temperature
!-- DENSX           dry air density (kg/m^3)
!-- PSTAR           P*=Psfc-Ptop
!-- QVS             water vapor mixing ratio (Kg/Kg)
!-- QCS             cloud mixing ratio (Kg/Kg)
!-- QIS             ice mixing ratio (Kg/Kg)
!-- DSIGFI          inverse of sigma layer delta
!-- G               gravity
!-- RD              gas constant for dry air (j/kg/k)
!-- CPAIR           specific heat of moist air (M^2 S^-2 K^-1)
!-- EDDYZ           eddy diffusivity KZ
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

      IMPLICIT NONE

!.......Arguments
  
!... Integer
      INTEGER,  INTENT(IN)   ::    its,ite, kts,kte,ims,ime,kms,kme
!... Real
      REAL , DIMENSION( ims:ime ),          INTENT(IN)  :: PBL, UST
      REAL ,                                INTENT(IN)  :: DTPBL, G, RD
      REAL , DIMENSION( kts:kte ),          INTENT(IN)  :: DSIGFI
      REAL , DIMENSION( its:ite ),          INTENT(IN)  :: MOL, PSTAR, CPAIR

      REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, TT,   &
                                                           QVS, QCS, QIS, DENSX
      REAL, DIMENSION( its:ite, kts:kte ), INTENT(IN) :: ZA, THETAV
      REAL, DIMENSION( its:ite, 0:kte )  , INTENT(IN) :: ZF
      
      REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: EDDYZ

!.......Local variables

!... Integer
      INTEGER  :: ILX, KL, KLM, K, I

!... Real
      REAL     :: ZOVL, PHIH, WT, ZSOL, ZFUNC, DZF, SS, GOTH, EDYZ
      REAL     :: RI, QMEAN, TMEAN, XLV, ALPH, CHI, ZK, SQL, DENSF, KZO
      REAL     :: FH
!... Parameters
      REAL, PARAMETER :: RV     = 461.5
      REAL, PARAMETER :: RC     = 0.25
      REAL, PARAMETER :: RLAM   = 80.0
      REAL, PARAMETER :: GAMH   = 16.0 !15.0  !  Holtslag and Boville (1993)
      REAL, PARAMETER :: BETAH  = 5.0   !  Holtslag and Boville (1993)
      REAL, PARAMETER :: KARMAN = 0.4
      REAL, PARAMETER :: EDYZ0  = 0.01  ! New Min Kz
!      REAL, PARAMETER :: EDYZ0  = 0.1
!--   IMVDIF      imvdif=1 for moist adiabat vertical diffusion
      INTEGER, PARAMETER :: imvdif = 1
!
      ILX = ite 
      KL  = kte
      KLM = kte - 1
      
      DO K = kts,KLM
        DO I = its,ILX
          EDYZ = 0.0
          ZOVL = 0.0
          DZF  = ZA(I,K+1) - ZA(I,K)
          KZO = EDYZ0
!--------------------------------------------------------------------------
          IF (ZF(I,K) .LT. PBL(I)) THEN
            ZOVL = ZF(I,K) / MOL(I)
            IF (ZOVL .LT. 0.0) THEN
              IF (ZF(I,K) .LT. 0.1 * PBL(I)) THEN
                PHIH = 1.0 / SQRT(1.0 - GAMH * ZOVL)
                WT   = UST(I) / PHIH
              ELSE
                ZSOL = 0.1 * PBL(I) / MOL(I)
                PHIH = 1.0 / SQRT(1.0 - GAMH * ZSOL)
                WT   = UST(I) / PHIH
              ENDIF
            ELSE IF (ZOVL .LT. 1.0) THEN
              PHIH = 1.0 + BETAH * ZOVL
              WT   = UST(I) / PHIH
            ELSE
              PHIH = BETAH + ZOVL
              WT   = UST(I) / PHIH
            ENDIF
            ZFUNC      = ZF(I,K) * (1.0 - ZF(I,K) / PBL(I)) ** 2
            EDYZ = KARMAN * WT * ZFUNC
          ENDIF
!--------------------------------------------------------------------------
          SS   = ((US(I,K+1) - US(I,K)) ** 2 + (VS(I,K+1) - VS(I,K)) ** 2)   &
                  / (DZF * DZF) + 1.0E-9
          GOTH = 2.0 * G / (THETAV(I,K+1) + THETAV(I,K))
          RI   = GOTH * (THETAV(I,K+1) - THETAV(I,K)) / (DZF * SS)
!--------------------------------------------------------------------------
!         Adjustment to vert diff in Moist air
          IF(imvdif.eq.1)then
            IF ((QCS(I,K)+QIS(I,K)) .GT. 0.01E-3 .OR. (QCS(I,K+1)+             &
                 QIS(I,K+1)) .GT. 0.01E-3) THEN
              QMEAN = 0.5 * (QVS(I,K) + QVS(I,K+1))
              TMEAN = 0.5 * (TT(I,K) + TT(I,K+1))
              XLV   = (2.501 - 0.00237 * (TMEAN - 273.15)) * 1.E6
              ALPH  =  XLV * QMEAN / RD / TMEAN
              CHI   =  XLV * XLV * QMEAN / CPAIR(I) / RV / TMEAN / TMEAN
              RI    = (1.0 + ALPH) * (RI -G * G / SS / TMEAN / CPAIR(I) *       &
                      ((CHI - ALPH) / (1.0 + CHI)))
            ENDIF
          ENDIF
!--------------------------------------------------------------------------
            
	        ZK  = 0.4 * ZF(I,K)
          SQL = (ZK * RLAM / (RLAM + ZK)) ** 2
            
          IF (RI .GE. 0.0) THEN
	          IF (ZF(I,K).LT.PBL(I).AND.ZOVL.GT.0.0) THEN
	            FH = MAX((1.-ZF(I,K)/PBL(I))**2,0.01) * PHIH **(-2)
                   SQL = ZK ** 2
	          ELSE
	            FH = (MAX(1.-RI/RC,0.01))**2
	          ENDIF
            EDDYZ(I,K) = KZO + SQRT(SS) * FH * SQL
          ELSE
            EDDYZ(I,K) = KZO + SQRT(SS * (1.0 - 25.0 * RI)) * SQL
          ENDIF
	  
          IF(EDYZ.GT.EDDYZ(I,K)) THEN
            EDDYZ(I,K) = EDYZ
          ENDIF

          EDDYZ(I,K) = MIN(1000.0,EDDYZ(I,K))
          EDDYZ(I,K) = MAX(KZO,EDDYZ(I,K))

          DENSF     = 0.5 * (DENSX(I,K+1) + DENSX(I,K))

          EDDYZ(I,K) = EDDYZ(I,K) * (DENSF * G / PSTAR(I)) ** 2 *       &
                       DTPBL * DSIGFI(K)*1E-6

        ENDDO             ! for I loop
      ENDDO               ! for k loop
!
      DO I = its,ILX
        EDDYZ(I,KL) = 0.0 ! EDDYZ(I,KLM) -- changed jp 3/08
      ENDDO

   END SUBROUTINE EDDYX
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!-------------------------------------------------------------------          
   SUBROUTINE ACM (DTPBL, PSTAR,  NOCONV, SIGMAF, DSIGH, DSIGHI, JX, &
                   KLPBL, PBL,   PBLSIG, MOL,  UST,                  &
                   TST, QST,  USTM,   EDDYZ,  DENSX,               &
                   US,    VS,     THETA,  QVS,    QCS,    QIS,     &
                   UX,    VX,     THETAX, QVX,    QCX,    QIX,     &
                   ids,ide, jds,jde, kds,kde,                      &
                   ims,ime, jms,jme, kms,kme,                      &
                   its,ite, jts,jte, kts,kte)
!**********************************************************************
!   PBL model called the Asymmetric Convective Model, Version 2 (ACM2) 
!   -- See top of module for summary and references
!
!---- REVISION HISTORY:
!   AX     3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM
!   JP and RG 8/2006 - updates
!
!**********************************************************************
!  ARGUMENTS:
!-- DTPBL           PBL time step
!-- PSTAR           Psurf - Ptop in cb
!-- NOCONV          If free convection =0, no; =1, yes
!-- SIGMAF          Sigma for full layer
!-- DSIGH           Sigma thickness
!-- DSIGHI          Inverse of sigma thickness
!-- JX              N-S index
!-- KLPBL           PBL level at K index
!-- PBL             PBL height in m
!-- PBLSIG          Sigma level for PBL 
!-- MOL             Monin-Obukhov length in 1D form
!-- UST             U* in 1D form
!-- TST             Theta* in 1D form
!-- QST             Q* in 1D form
!-- USTM            U* for computation of momemtum flux 
!-- EDDYZ           eddy diffusivity KZ
!-- DENSX           dry air density (kg/m^3)
!-- US              U wind 
!-- VS              V wind
!-- THETA           potential temperature
!-- QVS             water vapor mixing ratio (Kg/Kg)
!-- QCS             cloud mixing ratio (Kg/Kg)
!-- QIS             ice mixing ratio (Kg/Kg)
!-- UX              new U wind 
!-- VX              new V wind
!-- THETAX          new potential temperature
!-- QVX             new water vapor mixing ratio (Kg/Kg)
!-- QCX             new cloud mixing ratio (Kg/Kg)
!-- QIX             new ice mixing ratio (Kg/Kg)
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

      IMPLICIT NONE

!.......Arguments

!... Integer
      INTEGER,  INTENT(IN)      ::      ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, JX
      INTEGER,  DIMENSION( its:ite ), INTENT(IN)  :: NOCONV
      INTEGER,  DIMENSION( ims:ime ), INTENT(IN)  :: KLPBL

!... Real
      REAL , DIMENSION( ims:ime ),          INTENT(IN)  :: PBL, UST
      REAL ,                                INTENT(IN)  :: DTPBL
      REAL , DIMENSION( its:ite ),          INTENT(IN)  :: PSTAR, PBLSIG,  &
                                                           MOL, TST, &
                                                           QST, USTM
      REAL , DIMENSION( kts:kte ),          INTENT(IN)  :: DSIGHI, DSIGH
      REAL , DIMENSION( 0:kte ),            INTENT(IN)  :: SIGMAF
      REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT)  :: EDDYZ
      REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, THETA,   &
                                                           QVS, QCS, QIS, DENSX
      REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: UX, VX, THETAX,      &
                                                           QVX, QCX, QIX
!.......Local variables

!... Parameters
      INTEGER, PARAMETER :: NSP   = 6
!
!......ACM2 Parameters
!     INTEGER, PARAMETER :: IFACM = 0
!
      REAL,    PARAMETER :: G1000 = 9.8 * 1.0E-3
      REAL,    PARAMETER :: XX    = 0.5          ! FACTOR APPLIED TO CONV MIXING TIME STEP
      REAL,    PARAMETER :: KARMAN = 0.4

!... Integer
      INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L
      INTEGER :: KCBLMX
      INTEGER, DIMENSION( its:ite ) :: KCBL

!... Real
      REAL                               :: G1000I, MBMAX, HOVL, MEDDY, MBAR
      REAL                               :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1
      REAL, DIMENSION( its:ite )         :: PSTARI, FSACM, DTLIM
      REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN
      REAL, DIMENSION( 1:NSP, its:ite )  :: FS, BCBOTN
      REAL, DIMENSION( kts:kte )         :: XPLUS, XMINUS
      REAL  DELC
      REAL, DIMENSION( 1:NSP,its:ite,kts:kte  ) :: VCI

      REAL, DIMENSION( kts:kte )               :: AI, BI, CI, EI !, Y
      REAL, DIMENSION( 1:NSP,kts:kte )         :: DI, UI    
!
!--Start Exicutable ----

      ILX = ite
      KL  = kte
      KLM = kte - 1

      G1000I = 1.0 / G1000
      KCBLMX = 0
      MBMAX  = 0.0

!---COMPUTE ACM MIXING RATE
      DO I = its, ILX
        DTLIM(I)  = DTPBL
        PSTARI(I) = 1.0 / PSTAR(I)
        KCBL(I)   = 1
        FSACM(I)  = 0.0

        IF (NOCONV(I) .EQ. 1) THEN
          KCBL(I) = KLPBL(I)

!-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE
!--New couple ACM & EDDY-------------------------------------------------------------
          HOVL     = -PBL(I) / MOL(I)
          FSACM(I) = 1./(1.+((KARMAN/(HOVL))**0.3333)/(0.72*KARMAN))
          MEDDY    = EDDYZ(I,1) / (DTPBL * (PBLSIG(I) - SIGMAF(1)))
          MBAR     = MEDDY * FSACM(I)
          DO K = kts,KCBL(I)-1
            EDDYZ(I,K) = EDDYZ(I,K) * (1.0 - FSACM(I))
          ENDDO

          MBMAX = AMAX1(MBMAX,MBAR)
          DO K = kts+1,KCBL(I)
            MBARKS(K,I) = MBAR
            MDWN(K,I)   = MBAR * (PBLSIG(I) - SIGMAF(K-1)) * DSIGHI(K)
          ENDDO
          MBARKS(1,I) = MBAR
          MBARKS(KCBL(I),I) = MDWN(KCBL(I),I)
          MDWN(KCBL(I)+1,I) = 0.0
        ENDIF
      ENDDO                              ! end of I loop

      DO K = kts,KLM
        DO I = its,ILX
          EKZ   = EDDYZ(I,K) / DTPBL * DSIGHI(K)
          DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I))
        ENDDO
      ENDDO
       
      DO I = its,ILX 
        IF (NOCONV(I) .EQ. 1) THEN
          KCBLMX = AMAX0(KLPBL(I),KCBLMX)
          RZ     = (SIGMAF(KCBL(I)) - SIGMAF(1)) * DSIGHI(1)
          DTLIM(I)  = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I))
        ENDIF
      ENDDO

      DO K = kts,KL
        DO I = its,ILX
          VCI(1,I,K) = THETA(I,K)
          VCI(2,I,K) = QVS(I,K)
          VCI(3,I,K) = US(I,K)
          VCI(4,I,K) = VS(I,K)
          ! -- Also mix cloud water and ice IF necessary
          ! IF (IMOISTX.NE.1.AND.IMOISTX.NE.3) THEN  !!! Check other PBL models
          VCI(5,I,K) = QCS(I,K)
          VCI(6,I,K) = QIS(I,K)
        ENDDO
      ENDDO

      NSPX=6

      DO I = its,ILX
        FS(1,I) = -UST(I) * TST(I) * DENSX(I,1) * PSTARI(I)
        FS(2,I) = -UST(I) * QST(I) * DENSX(I,1) * PSTARI(I)
        FM      = -USTM(I) * USTM(I) * DENSX(I,1) * PSTARI(I)
        WSPD    = SQRT(US(I,1) * US(I,1) + VS(I,1) * VS(I,1)) + 1.E-9
        FS(3,I) = FM * US(I,1) / WSPD
        FS(4,I) = FM * VS(I,1) / WSPD
        FS(5,I) = 0.0
        FS(6,I) = 0.0                      ! SURFACE FLUXES OF CLOUD WATER AND ICE = 0
      ENDDO
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      DO I = its,ILX      

        NLP   = INT(DTPBL / DTLIM(I) + 1.0)
        DTS   = (DTPBL / NLP)
        DTRAT = DTS / DTPBL
        DO NL = 1,NLP           ! LOOP OVER SUB TIME LOOP              

!-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES

          DO K = kts,KL
            AI(K) = 0.0
            BI(K) = 0.0
            CI(K) = 0.0
            EI(K) = 0.0
          ENDDO

          DO K = 2, KCBL(I)
            EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DSIGH(K) * DSIGHI(K-1)
            BI(K)   = 1.0 + CRANKP * MDWN(K,I) * DTS
            AI(K)   = -CRANKP * MBARKS(K,I) * DTS
          ENDDO

          EI(1) = EI(1) -EDDYZ(I,1) * CRANKP * DSIGHI(1 )* DTRAT
          AI(2) = AI(2) -EDDYZ(I,1) * CRANKP * DSIGHI(2) * DTRAT

          DO K =  KCBL(I)+1, KL
            BI(K) = 1.0
          ENDDO

          DO K = 2,KL
            XPLUS(K)  = EDDYZ(I,K) * DSIGHI(K) * DTRAT
            XMINUS(K) = EDDYZ(I,K-1) * DSIGHI(K) * DTRAT
            CI(K)     = - XMINUS(K) * CRANKP
            EI(K)     = EI(K) - XPLUS(K) * CRANKP
            BI(K)     = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP
          ENDDO

          IF (NOCONV(I) .EQ. 1) THEN
            BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBLSIG(I) - SIGMAF(1)) *    &
                    DTS * DSIGHI(1) + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
          ELSE
            BI(1) = 1.0  + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
          ENDIF


          DO K = 1,KL
            DO L = 1,NSPX                    
              DI(L,K) = 0.0
            ENDDO
          ENDDO
!
!**   COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION
          DO K = 2,KCBL(I)
            DO L = 1,NSPX                    
              DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) *          &
                 VCI(L,I,K) + DSIGH(K+1) * DSIGHI(K) *                  &
                        MDWN(K+1,I) * VCI(L,I,K+1))
              DI(L,K)   = VCI(L,I,K) + (1.0 - CRANKP) * DELC
            ENDDO
          ENDDO

          DO K = KCBL(I)+1, KL
            DO L = 1,NSPX                    
              DI(L,K) = VCI(L,I,K)
            ENDDO
          ENDDO

          DO K = 2,KL
            IF (K .EQ. KL) THEN
              DO L = 1,NSPX                    
                DI(L,K) = DI(L,K)  - (1.0 - CRANKP) * XMINUS(K) *                  &
                          (VCI(L,I,K) - VCI(L,I,K-1))
              ENDDO
            ELSE
              DO L = 1,NSPX                    
                DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) *                   &
                          (VCI(L,I,K+1) - VCI(L,I,K))  -                         &
                          (1.0 - CRANKP) * XMINUS(K) *                           &
                          (VCI(L,I,K) - VCI(L,I,K-1))
              ENDDO
            ENDIF
          ENDDO

          IF (NOCONV(I) .EQ. 1) THEN
            DO L = 1,NSPX                    
              F1    = -G1000I * (MBARKS(1,I) *                                &
                      (PBLSIG(I) - SIGMAF(1)) * VCI(L,I,1) -                  &
                      MDWN(2,I) * VCI(L,I,2) * DSIGH(2))

              DI(L,1) = VCI(L,I,1) - G1000 * (FS(L,I) - (1.0 - CRANKP)        &
                        * F1) * DSIGHI(1) * DTS
            ENDDO
          ELSE
            DO L = 1,NSPX                    
              DI(L,1) = VCI(L,I,1) - G1000 * FS(L,I) * DSIGHI(1) * DTS
            ENDDO
          ENDIF
          DO L = 1,NSPX                    
            DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZ(I,1) * DSIGHI(1)      &
                     * DTRAT * (VCI(L,I,2) - VCI(L,I,1))
          ENDDO
          IF ( NOCONV(I) .EQ. 1 ) THEN
            CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX)
          ELSE
            CALL TRI (CI, BI, EI, DI, UI, KL, NSPX)
          END IF
!
!-- COMPUTE NEW THETAV AND Q
          DO K = 1,KL
            DO L = 1,NSPX                    
              VCI(L,I,K) = UI(L,K)
            ENDDO
          ENDDO

        ENDDO                   ! END I LOOP
      ENDDO                     ! END SUB TIME LOOP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
      DO K = kts,KL
        DO I = its,ILX
          THETAX(I,K) = VCI(1,I,K)
          QVX(I,K)    = VCI(2,I,K)
          UX(I,K)     = VCI(3,I,K)
          VX(I,K)     = VCI(4,I,K)
        ENDDO
      ENDDO

      DO K = kts,KL
        DO I = its,ILX
          QCX(I,K) = VCI(5,I,K)
          QIX(I,K) = VCI(6,I,K)
        ENDDO
      ENDDO

   END SUBROUTINE ACM
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   SUBROUTINE MATRIX(A,B,C,D,E,X,KL,NSP)
   
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   IMPLICIT NONE
!
!-- Bordered band diagonal matrix solver for ACM2

!-- ACM2 Matrix is in this form:
!   B1 E1
!   A2 B2 E2
!   A3 C3 B3 E3
!   A4    C4 B4 E4
!   A5       C5 B5 E5
!   A6          C6 B6

!--Upper Matrix is
!  U11 U12
!      U22 U23
!          U33 U34
!              U44 U45
!                  U55 U56
!                      U66

!--Lower Matrix is:
!  1
! L21  1
! L31 L32  1
! L41 L42 L43  1
! L51 L52 L53 L54  1
! L61 L62 L63 L64 L65 1
!---------------------------------------------------------
!...Arguments
      INTEGER, INTENT(IN)   :: KL
      INTEGER, INTENT(IN)   :: NSP
      REAL A(KL),B(KL),E(KL)
      REAL C(KL),D(NSP,KL),X(NSP,KL)

!...Locals
      REAL Y(NSP,KL),AIJ,SUM
      REAL L(KL,KL),UII(KL),UIIP1(KL),RUII(KL)
      INTEGER I,J,V

!-- Define Upper and Lower matrices
      L(1,1) = 1.
      UII(1) = B(1)
      RUII(1) = 1./UII(1)
      DO I = 2, KL
	      L(I,I) = 1.
	      L(I,1) = A(I)/B(1)
        UIIP1(I-1)=E(I-1)
	      IF(I.GE.3) THEN
	        DO J = 2,I-1
	          IF(I.EQ.J+1) THEN
	            AIJ = C(I)
	          ELSE
	            AIJ = 0.
	          ENDIF
	          L(I,J) = (AIJ-L(I,J-1)*E(J-1))/      &
                      (B(J)-L(J,J-1)*E(J-1))
	        ENDDO
	      ENDIF
      ENDDO
	  
      DO I = 2,KL
        UII(I) = B(I)-L(I,I-1)*E(I-1)
        RUII(I) = 1./UII(I)
      ENDDO
  
!-- Forward sub for Ly=d
      DO V= 1, NSP
        Y(V,1) = D(V,1)
        DO I=2,KL
	        SUM = D(V,I)
	        DO J=1,I-1
	          SUM = SUM - L(I,J)*Y(V,J)
	        ENDDO
	        Y(V,I) = SUM
        ENDDO
      ENDDO

!-- Back sub for Ux=y

      DO V= 1, NSP
        X(V,KL) = Y(V,KL)*RUII(KL)
      ENDDO
      DO I = KL-1,1,-1
        DO V= 1, NSP
         X(V,I) = (Y(V,I)-UIIP1(I)*X(V,I+1))*RUII(I)
        ENDDO
      ENDDO

   END SUBROUTINE MATRIX


!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
      SUBROUTINE TRI ( L, D, U, B, X,KL,NSP)
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

!  FUNCTION:
!    Solves tridiagonal system by Thomas algorithm. 
!   The associated tri-diagonal system is stored in 3 arrays
!   D : diagonal
!   L : sub-diagonal
!   U : super-diagonal
!   B : right hand side function
!   X : return solution from tridiagonal solver

!     [ D(1) U(1) 0    0    0 ...       0     ]
!     [ L(2) D(2) U(2) 0    0 ...       .     ]
!     [ 0    L(3) D(3) U(3) 0 ...       .     ]
!     [ .       .     .     .           .     ] X(i) = B(i)
!     [ .             .     .     .     0     ]
!     [ .                   .     .     .     ]
!     [ 0                           L(n) D(n) ]

!-----------------------------------------------------------------------

      IMPLICIT NONE

! Arguments:

      INTEGER, INTENT(IN)   :: KL
      INTEGER, INTENT(IN)   :: NSP

      REAL        L( KL )               ! subdiagonal
      REAL        D(KL)   ! diagonal
      REAL        U( KL )               ! superdiagonal
      REAL        B(NSP,KL )   ! R.H. side
      REAL        X( NSP,KL )   ! solution

! Local Variables:

      REAL        GAM( KL )
      REAL        BET
      INTEGER     V, K

! Decomposition and forward substitution:
      BET = 1.0 / D( 1 )
      DO V = 1, NSP
         X( V,1 ) = BET * B(V,1 )
      ENDDO

      DO K = 2, KL
        GAM(K ) = BET * U( K-1 )
        BET = 1.0 / ( D( K ) - L( K ) * GAM( K ) )
	      DO V = 1, NSP
           X( V, K ) = BET * ( B( V,K ) - L( K ) * X( V,K-1 ) )
	      ENDDO
      ENDDO

! Back-substitution:

      DO K = KL - 1, 1, -1
        DO V = 1, NSP
          X( V,K ) = X( V,K ) - GAM( K+1 ) * X( V,K+1 )
        ENDDO
      ENDDO
     
  END SUBROUTINE TRI
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

END MODULE module_bl_acm
                        
Back to Top