phys/module_bl_mynn_wrapper.F

   1 !> \file module_bl_mynn_wrapper.F90
   2 !!  This serves as the interface between the WRF PBL driver and the MYNN
   3 !!  eddy-diffusivity mass-flux scheme in module_bl_mynn.F.
   4
   5 !>\ingroup gsd_mynn_edmf
   6 !> The following references best describe the code within
   7 !!    Olson et al. (2019, NOAA Technical Memorandum)
   8 !!    Nakanishi and Niino (2009) \cite NAKANISHI_2009
   9       MODULE module_bl_mynn_wrapper
  10
  11       use module_bl_mynn_common
  12
  13       contains
  14
  15 !> \section arg_table_mynnedmf_wrapper_init Argument Table
  16 !! \htmlinclude mynnedmf_wrapper_init.html
  17 !!
  18       subroutine mynnedmf_wrapper_init (              &
  19         &  RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN,&
  20         &  RQIBLTEN,QKE,                              &
  21         &  restart,allowed_to_read,                   &
  22         &  P_QC,P_QI,PARAM_FIRST_SCALAR,              &
  23         &  IDS,IDE,JDS,JDE,KDS,KDE,                   &
  24         &  IMS,IME,JMS,JME,KMS,KME,                   &
  25         &  ITS,ITE,JTS,JTE,KTS,KTE                    )
  26
  27         implicit none
  28
  29         LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
  30
  31         INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,  &
  32              &                IMS,IME,JMS,JME,KMS,KME,  &
  33              &                ITS,ITE,JTS,JTE,KTS,KTE
  34
  35
  36         REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: &
  37              &RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,                 &
  38              &RQCBLTEN,RQIBLTEN,QKE
  39
  40         INTEGER,  intent(in) :: P_QC,P_QI,PARAM_FIRST_SCALAR
  41
  42         INTEGER :: I,J,K,ITF,JTF,KTF
  43
  44         JTF=MIN0(JTE,JDE-1)
  45         KTF=MIN0(KTE,KDE-1)
  46         ITF=MIN0(ITE,IDE-1)
  47
  48         IF (.NOT.RESTART) THEN
  49          DO J=JTS,JTF
  50           DO K=KTS,KTF
  51            DO I=ITS,ITF
  52               RUBLTEN(i,k,j)=0.
  53               RVBLTEN(i,k,j)=0.
  54               RTHBLTEN(i,k,j)=0.
  55               RQVBLTEN(i,k,j)=0.
  56               if( p_qc >= param_first_scalar ) RQCBLTEN(i,k,j)=0.
  57               if( p_qi >= param_first_scalar ) RQIBLTEN(i,k,j)=0.
  58            ENDDO
  59           ENDDO
  60          ENDDO
  61         ENDIF
  62
  63       end subroutine mynnedmf_wrapper_init
  64
  65       subroutine mynnedmf_wrapper_finalize ()
  66       end subroutine mynnedmf_wrapper_finalize
  67
  68 ! \brief This scheme (1) performs pre-mynnedmf work, (2) runs the mynnedmf, and (3) performs post-mynnedmf work
  69 !> \section arg_table_mynnedmf_wrapper_run Argument Table
  70 !! \htmlinclude mynnedmf_wrapper_run.html
  71 !!
  72 SUBROUTINE mynnedmf_wrapper_run(                 &
  73      &  initflag,restart,cycling,                &
  74      &  delt,dz,dxc,znt,                         &
  75      &  u,v,w,th,                                &
  76      &  qv,qc,qi,qnc,qni,qnwfa,qnifa,qnbca,      &
  77 !     &  ozone,                                   &
  78      &  p,exner,rho,t3d,                         &
  79      &  xland,ts,qsfc,ps,                        &
  80      &  ust,ch,hfx,qfx,rmol,wspd,                &
  81      &  uoce,voce,                               &
  82      &  qke,qke_adv,sh3d,sm3d,                   &
  83 !--- chem/smoke
  84 #if (WRF_CHEM == 1)
  85      &  mix_chem,chem3d,vd3d,nchem,kdvel,        &
  86      &  ndvel,num_vert_mix,                      &
  87 !     &  frp_mean,emis_ant_no,enh_mix,            & !to be included soon
  88 #endif
  89 !--- end chem/smoke
  90      &  Tsq,Qsq,Cov,                             &
  91      &  rublten,rvblten,rthblten,                &
  92      &  rqvblten,rqcblten,rqiblten,              &
  93      &  rqncblten,rqniblten,                     &
  94      &  rqnwfablten,rqnifablten,rqnbcablten,     &
  95 !     &  ro3blten,                                &
  96      &  exch_h,exch_m,pblh,kpbl,el_pbl,          &
  97      &  dqke,qwt,qshear,qbuoy,qdiss,             &
  98      &  qc_bl,qi_bl,cldfra_bl,                   &
  99      &  edmf_a,edmf_w,edmf_qt,                   &
 100      &  edmf_thl,edmf_ent,edmf_qc,               &
 101      &  sub_thl3d,sub_sqv3d,                     &
 102      &  det_thl3d,det_sqv3d,                     &
 103      &  nupdraft,maxMF,ktop_plume,               &
 104      &  rthraten,                                &
 105      &  tke_budget,         bl_mynn_tkeadvect,   &
 106      &  bl_mynn_cloudpdf,   bl_mynn_mixlength,   &
 107      &  icloud_bl,          bl_mynn_edmf,        &
 108      &  bl_mynn_edmf_mom,   bl_mynn_edmf_tke,    &
 109      &  bl_mynn_cloudmix,   bl_mynn_mixqt,       &
 110      &  bl_mynn_output,     bl_mynn_closure,     &
 111      &  bl_mynn_mixscalars,                      &
 112      &  spp_pbl,pattern_spp_pbl,                 &
 113      &  flag_qc,flag_qi,                         &
 114      &  flag_qnc,flag_qni,                       &
 115      &  flag_qnwfa,flag_qnifa,flag_qnbca,        &
 116      &  ids,ide,jds,jde,kds,kde,                 &
 117      &  ims,ime,jms,jme,kms,kme,                 &
 118      &  its,ite,jts,jte,kts,kte                  )
 119
 120 !     use module_gfs_machine,  only : kind_phys
 121      use module_bl_mynn, only: mynn_bl_driver
 122
 123 !-------------------------------------------------------------------
 124      implicit none
 125 !-------------------------------------------------------------------
 126
 127      !smoke/chem: disclaimer: all smoke-related variables are still
 128      !considered under development in CCPP. Until that work is
 129      !completed, these flags/arrays must be kept hard-coded as is.
 130 #if (WRF_CHEM == 1)
 131      logical, intent(in) :: mix_chem
 132      integer, intent(in) :: nchem, ndvel, kdvel, num_vert_mix
 133      logical, parameter ::                                  &
 134      &       rrfs_sd    =.false.,                           &
 135      &       smoke_dbg  =.false.,                           &
 136      &       enh_mix    =.false.
 137 #else
 138      logical, parameter ::                                  &
 139      &       mix_chem   =.false.,                           &
 140      &       enh_mix    =.false.,                           &
 141      &       rrfs_sd    =.false.,                           &
 142      &       smoke_dbg  =.false.
 143      integer, parameter :: nchem=2, ndvel=2, kdvel=1,       &
 144              num_vert_mix = 1
 145 #endif
 146
 147 ! NAMELIST OPTIONS (INPUT):
 148      logical, intent(in) ::                                 &
 149      &       bl_mynn_tkeadvect,                             &
 150      &       cycling
 151       integer, intent(in) ::                                &
 152      &       bl_mynn_cloudpdf,                              &
 153      &       bl_mynn_mixlength,                             &
 154      &       icloud_bl,                                     &
 155      &       bl_mynn_edmf,                                  &
 156      &       bl_mynn_edmf_mom,                              &
 157      &       bl_mynn_edmf_tke,                              &
 158      &       bl_mynn_cloudmix,                              &
 159      &       bl_mynn_mixqt,                                 &
 160      &       bl_mynn_output,                                &
 161      &       bl_mynn_mixscalars,                            &
 162      &       spp_pbl,                                       &
 163      &       tke_budget
 164       real, intent(in) ::                                   &
 165      &       bl_mynn_closure
 166
 167       logical, intent(in) ::                                &
 168      &       FLAG_QI, FLAG_QNI, FLAG_QC, FLAG_QNC,          &
 169      &       FLAG_QNWFA, FLAG_QNIFA, FLAG_QNBCA
 170       logical, parameter :: FLAG_OZONE = .false.
 171
 172 !MYNN-1D
 173       REAL,    intent(in) :: delt, dxc
 174       LOGICAL, intent(in) :: restart
 175       INTEGER :: i, j, k, itf, jtf, ktf, n
 176       INTEGER, intent(in) :: initflag,                      &
 177      &           IDS,IDE,JDS,JDE,KDS,KDE,                   &
 178      &           IMS,IME,JMS,JME,KMS,KME,                   &
 179      &           ITS,ITE,JTS,JTE,KTS,KTE
 180
 181 !MYNN-3D
 182       real, dimension(ims:ime,kms:kme,jms:jme), intent(in) ::          &
 183      &       u,v,w,t3d,th,rho,exner,p,dz
 184       real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) ::       &
 185      &       rublten,rvblten,rthblten,                                 &
 186      &       rqvblten,rqcblten,rqiblten,                               &
 187      &       rqncblten,rqniblten,                                      &
 188      &       rqnwfablten,rqnifablten,rqnbcablten !,ro3blten
 189       real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) ::       &
 190      &       qke, qke_adv, el_pbl, sh3d, sm3d
 191       real, dimension(ims:ime,kms:kme,jms:jme), intent(inout) ::       &
 192      &       Tsq, Qsq, Cov, exch_h, exch_m
 193       real, dimension(ims:ime,kms:kme,jms:jme), intent(in) :: rthraten
 194
 195 !optional 3D arrays
 196       real, dimension(ims:ime,kms:kme,jms:jme), optional, intent(in)  ::    &
 197      &       pattern_spp_pbl
 198       real, dimension(ims:ime,kms:kme,jms:jme), optional, intent(inout) ::  &
 199      &       qc_bl, qi_bl, cldfra_bl
 200       real, dimension(ims:ime,kms:kme,jms:jme), optional, intent(inout) ::  &
 201      &       edmf_a,edmf_w,edmf_qt,                                    &
 202      &       edmf_thl,edmf_ent,edmf_qc,                                &
 203      &       sub_thl3d,sub_sqv3d,det_thl3d,det_sqv3d
 204       real, dimension(ims:ime,kms:kme,jms:jme), optional, intent(inout) ::  &
 205      &       dqke,qWT,qSHEAR,qBUOY,qDISS
 206       real, dimension(ims:ime,kms:kme,jms:jme), optional, intent(inout) ::  &
 207      &       qv,qc,qi,qnc,qni,qnwfa,qnifa,qnbca!,o3
 208
 209 !optional 2D arrays for passing into module_bl_myn.F
 210       real, allocatable, dimension(:,:) ::                             &
 211      &       qc_bl2d, qi_bl2d, cldfra_bl2d, pattern_spp_pbl2d
 212       real, allocatable, dimension(:,:) ::                             &
 213      &       edmf_a2d,edmf_w2d,edmf_qt2d,                              &
 214      &       edmf_thl2d,edmf_ent2d,edmf_qc2d,                          &
 215      &       sub_thl2d,sub_sqv2d,det_thl2d,det_sqv2d
 216       real, allocatable, dimension(:,:) ::                             &
 217      &       dqke2d,qWT2d,qSHEAR2d,qBUOY2d,qDISS2d
 218       real, allocatable, dimension(:,:) ::                             &
 219      &       qc2d,qi2d,qnc2d,qni2d,qnwfa2d,qnifa2d,qnbca2d!,o32d
 220
 221 !smoke/chem arrays - no if-defs in module_bl_mynn.F, so must define arrays
 222 #if (WRF_CHEM == 1)
 223       real, dimension(ims:ime,kms:kme,jms:jme,nchem), intent(in) :: chem3d
 224       real, dimension(ims:ime,kdvel,jms:jme, ndvel), intent(in)  :: vd3d
 225       real, dimension(ims:ime,kms:kme,nchem) :: chem
 226       real, dimension(ims:ime,ndvel)         :: vd
 227       real, dimension(ims:ime)               :: frp_mean, emis_ant_no
 228 #else
 229       real, dimension(ims:ime,kms:kme,nchem) :: chem
 230       real, dimension(ims:ime,ndvel)         :: vd
 231       real, dimension(ims:ime) :: frp_mean, emis_ant_no
 232 #endif
 233
 234 !MYNN-2D
 235       real, dimension(ims:ime,jms:jme), intent(in) ::                  &
 236      &        xland,ts,qsfc,ps,ch
 237       real, dimension(ims:ime,jms:jme), intent(inout) ::               &
 238      &        znt,pblh,maxmf,rmol,hfx,qfx,ust,wspd,                    &
 239      &        uoce,voce
 240       integer, dimension(ims:ime,jms:jme), intent(inout) ::            &
 241      &        kpbl,nupdraft,ktop_plume
 242
 243 !Local
 244       real, dimension(ims:ime,kms:kme)         :: delp,sqv,sqc,sqi
 245       real, dimension(ims:ime)                 :: dx
 246       logical, parameter                       :: debug = .false.
 247       real, dimension(ims:ime,kms:kme,jms:jme) :: ozone,r03blten
 248
 249    !write(0,*)"=============================================="
 250    !write(0,*)"in mynn wrapper..."
 251    !write(0,*)"initflag=",initflag
 252    !write(0,*)"restart =",restart
 253
 254    jtf=MIN0(JTE,JDE-1)
 255    ktf=MIN0(KTE,KDE-1)
 256    itf=MIN0(ITE,IDE-1)
 257
 258    !For now, initialized bogus array
 259    ozone=0.0
 260    r03blten=0.0
 261
 262    !Allocate any arrays being used
 263    if (icloud_bl > 0) then
 264       allocate(qc_bl2d(ims:ime,kms:kme))
 265       allocate(qi_bl2d(ims:ime,kms:kme))
 266       allocate(cldfra_bl2d(ims:ime,kms:kme))
 267       qc_bl2d=0.0
 268       qi_bl2d=0.0
 269       cldfra_bl2d=0.0
 270    endif
 271    if (spp_pbl > 0) then
 272       allocate(pattern_spp_pbl2d(ims:ime,kms:kme))
 273    endif
 274    if (bl_mynn_output > 0) then
 275       allocate(edmf_a2d(ims:ime,kms:kme))
 276       allocate(edmf_w2d(ims:ime,kms:kme))
 277       allocate(edmf_qt2d(ims:ime,kms:kme))
 278       allocate(edmf_thl2d(ims:ime,kms:kme))
 279       allocate(edmf_ent2d(ims:ime,kms:kme))
 280       allocate(edmf_qc2d(ims:ime,kms:kme))
 281       allocate(sub_thl2d(ims:ime,kms:kme))
 282       allocate(sub_sqv2d(ims:ime,kms:kme))
 283       allocate(det_thl2d(ims:ime,kms:kme))
 284       allocate(det_sqv2d(ims:ime,kms:kme))
 285    endif
 286    if (tke_budget .eq. 1) then
 287       allocate(dqke2d(ims:ime,kms:kme))
 288       allocate(qWT2d(ims:ime,kms:kme))
 289       allocate(qSHEAR2d(ims:ime,kms:kme))
 290       allocate(qBUOY2d(ims:ime,kms:kme))
 291       allocate(qDISS2d(ims:ime,kms:kme))
 292       dqke2d  =0.0
 293       qWT2d   =0.0
 294       qSHEAR2d=0.0
 295       qBUOY2d =0.0
 296       qDISS2d =0.0
 297    endif
 298    if (flag_qc) then
 299       allocate(qc2d(ims:ime,kms:kme))
 300       qc2d=0.0
 301    endif
 302    if (flag_qi) then
 303       allocate(qi2d(ims:ime,kms:kme))
 304       qi2d=0.0
 305    endif
 306    if (flag_qnc) then
 307       allocate(qnc2d(ims:ime,kms:kme))
 308       qnc2d=0.0
 309    endif
 310    if (flag_qni) then
 311       allocate(qni2d(ims:ime,kms:kme))
 312       qni2d=0.0
 313    endif
 314    if (flag_qnwfa) then
 315       allocate(qnwfa2d(ims:ime,kms:kme))
 316       qnwfa2d=0.0
 317    endif
 318    if (flag_qnifa) then
 319       allocate(qnifa2d(ims:ime,kms:kme))
 320       qnifa2d=0.0
 321    endif
 322    if (flag_qnbca) then
 323       allocate(qnbca2d(ims:ime,kms:kme))
 324       qnbca2d=0.0
 325    endif
 326    !---------------------------------
 327    !Begin looping in the j-direction
 328    !---------------------------------
 329    do j = jts,jtf
 330
 331       !need sgs cloud info input for diagnostic-decay
 332       if (icloud_bl > 0) then
 333          do k=kts,ktf
 334          do i=its,itf
 335             qc_bl2d(i,k)     = qc_bl(i,k,j)
 336             qi_bl2d(i,k)     = qi_bl(i,k,j)
 337             cldfra_bl2d(i,k) = cldfra_bl(i,k,j)
 338          enddo
 339          enddo
 340       endif
 341
 342       !spp input
 343       if (spp_pbl > 0) then
 344          do k=kts,ktf
 345          do i=its,itf
 346             pattern_spp_pbl2d(i,k) = pattern_spp_pbl(i,k,j)
 347          enddo
 348          enddo
 349       endif
 350
 351       !intialize moist species
 352       if (flag_qc) then
 353          do k=kts,ktf
 354          do i=its,itf
 355             qc2d(i,k) = qc(i,k,j)
 356          enddo
 357          enddo
 358       endif
 359       if (flag_qi) then
 360          do k=kts,ktf
 361          do i=its,itf
 362             qi2d(i,k) = qi(i,k,j)
 363          enddo
 364          enddo
 365       endif
 366       if (flag_qnc) then
 367          do k=kts,ktf
 368          do i=its,itf
 369             qnc2d(i,k) = qnc(i,k,j)
 370          enddo
 371          enddo
 372       endif
 373       if (flag_qni) then
 374          do k=kts,ktf
 375          do i=its,itf
 376             qni2d(i,k) = qni(i,k,j)
 377          enddo
 378          enddo
 379       endif
 380       if (flag_qnwfa) then
 381          do k=kts,ktf
 382          do i=its,itf
 383             qnwfa2d(i,k) = qnwfa(i,k,j)
 384          enddo
 385          enddo
 386       endif
 387       if (flag_qnifa) then
 388          do k=kts,ktf
 389          do i=its,itf
 390             qnifa2d(i,k) = qnifa(i,k,j)
 391          enddo
 392          enddo
 393       endif
 394       if (flag_qnbca) then
 395          do k=kts,ktf
 396          do i=its,itf
 397             qnbca2d(i,k) = qnbca(i,k,j)
 398          enddo
 399          enddo
 400       endif
 401
 402 #if (WRF_CHEM == 1)
 403       if (mix_chem) then
 404          do n=1,nchem
 405          do k=kts,ktf
 406          do i=its,itf
 407             chem(i,k,n)=chem3d(i,k,j,n)
 408          enddo
 409          enddo
 410          enddo
 411
 412          !set kdvel =1
 413          do n=1,ndvel
 414          do i=its,itf
 415             vd(i,n)  =vd3d(i,1,j,n)
 416          enddo
 417          enddo
 418       endif
 419       frp_mean = 0.0
 420       emis_ant_no = 0.0
 421 #else
 422       chem     = 0.0
 423       vd       = 0.0
 424       frp_mean = 0.0
 425       emis_ant_no = 0.0
 426 #endif
 427
 428       ! Check incoming moist species to ensure non-negative values
 429       ! First, create pressure differences (delp) across model layers
 430       do i=its,itf
 431          dx(i)=dxc
 432 !         delp(i,1)  = ps(i,j) - (p(i,2,j)*dz(i,1,j) + p(i,1,j)*dz(i,2,j))/(dz(i,1,j)+dz(i,2,j))
 433 !         do k=2,kte-1
 434 !            delp(i,k) = (p(i,k,j)*dz(i,k-1,j) + p(i,k-1,j)*dz(i,k,j))/(dz(i,k,j)+dz(i,k-1,j)) - &
 435 !                        (p(i,k+1,j)*dz(i,k,j) + p(i,k,j)*dz(i,k+1,j))/(dz(i,k,j)+dz(i,k+1,j))
 436 !         enddo
 437 !         delp(i,kte) = delp(i,kte-1)
 438       enddo
 439
 440 !      do i=its,itf
 441 !         call moisture_check2(kte, delt,                 &
 442 !                              delp(i,:), exner(i,:,j),   &
 443 !                              qv(i,:,j), qc(i,:,j),      &
 444 !                              qi(i,:,j), t3d(i,:,j)      )
 445 !      enddo
 446
 447       !In WRF, mixing ratio is incoming. Convert to specific humidity:
 448        do k=kts,ktf
 449           do i=its,itf
 450              sqv(i,k)=qv(i,k,j)/(1.0 + qv(i,k,j))
 451              sqc(i,k)=qc2d(i,k)/(1.0 + qv(i,k,j))
 452              sqi(i,k)=qi2d(i,k)/(1.0 + qv(i,k,j))
 453          enddo
 454       enddo
 455
 456 !      do i=its,ite
 457 !         ts(i,j)=tsurf(i,j)/exner(i,1,j)  !theta
 458 !      enddo
 459
 460       if (debug) then
 461          print*
 462          write(0,*)"===CALLING mynn_bl_driver; input:"
 463          print*,"tke_budget=",tke_budget
 464          print*,"bl_mynn_tkeadvect=",bl_mynn_tkeadvect
 465          print*,"bl_mynn_cloudpdf=",bl_mynn_cloudpdf
 466          print*,"bl_mynn_mixlength=",bl_mynn_mixlength
 467          print*,"bl_mynn_edmf=",bl_mynn_edmf
 468          print*,"bl_mynn_edmf_mom=",bl_mynn_edmf_mom
 469          print*,"bl_mynn_edmf_tke=",bl_mynn_edmf_tke
 470          print*,"bl_mynn_cloudmix=",bl_mynn_cloudmix
 471          print*,"bl_mynn_mixqt=",bl_mynn_mixqt
 472          print*,"icloud_bl=",icloud_bl
 473          print*,"T:",t3d(its,1,j),t3d(its,2,j),t3d(its,kte,j)
 474          print*,"TH:",th(its,1,j),th(its,2,j),th(its,kte,j)
 475          print*,"rho:",rho(its,1,j),rho(its,2,j),rho(its,kte,j)
 476          print*,"exner:",exner(its,1,j),exner(its,2,j),exner(its,kte,j)
 477          print*,"p:",p(its,1,j),p(its,2,j),p(its,kte,j)
 478          print*,"dz:",dz(its,1,j),dz(its,2,j),dz(its,kte,j)
 479          print*,"u:",u(its,1,j),u(its,2,j),u(its,kte,j)
 480          print*,"v:",v(its,1,j),v(its,2,j),v(its,kte,j)
 481          print*,"sqv:",sqv(its,1),sqv(its,2),sqv(its,kte)
 482          print*,"sqc:",sqc(its,1),sqc(its,2),sqc(its,kte)
 483          print*,"sqi:",sqi(its,1),sqi(its,2),sqi(its,kte)
 484          print*,"rmol:",rmol(its,j)," ust:",ust(its,j)
 485          print*,"dx=",dx(its),"initflag=",initflag
 486          print*,"Thetasurf:",ts(its,j)
 487          print*,"HFX:",hfx(its,j)," qfx",qfx(its,j)
 488          print*,"qsfc:",qsfc(its,j)," ps:",ps(its,j)
 489          print*,"wspd:",wspd(its,j)
 490          print*,"znt:",znt(its,j)," delt=",delt
 491          print*,"ite=",ite," kte=",kte
 492          print*,"PBLH=",pblh(its,j)," KPBL=",KPBL(its,j)," xland=",xland(its,j)
 493          print*," ch=",ch(its,j)
 494          print*,"qke:",qke(its,1,j),qke(its,2,j),qke(its,kte,j)
 495          print*,"el_pbl:",el_pbl(its,1,j),el_pbl(its,2,j),el_pbl(its,kte,j)
 496          print*,"Sh3d:",Sh3d(its,1,j),sh3d(its,2,j),sh3d(its,kte,j)
 497          print*,"max cf_bl:",maxval(cldfra_bl(its,:,j))
 498       endif
 499
 500 !print*,"In mynn wrapper, calling mynn_bl_driver"
 501               CALL  mynn_bl_driver(                                    &
 502      &        initflag=initflag,restart=restart,cycling=cycling,       &
 503      &        delt=delt,dz=dz(:,:,j),dx=dx,znt=znt(:,j),               &
 504      &        u=u(:,:,j),v=v(:,:,j),w=w(:,:,j),                        &
 505      &        th=th(:,:,j),sqv3D=sqv,sqc3D=sqc,                        &
 506      &        sqi3D=sqi,qnc=qnc2d,qni=qni2d,                           &
 507      &        qnwfa=qnwfa2d,qnifa=qnifa2d,                             &
 508      &        ozone=ozone(:,:,j),                                      &
 509      &        p=p(:,:,j),exner=exner(:,:,j),rho=rho(:,:,j),            &
 510      &        T3D=t3d(:,:,j),xland=xland(:,j),                         &
 511      &        ts=ts(:,j),qsfc=qsfc(:,j),ps=ps(:,j),                    &
 512      &        ust=ust(:,j),ch=ch(:,j),hfx=hfx(:,j),qfx=qfx(:,j),       &
 513      &        rmol=rmol(:,j),wspd=wspd(:,j),                           &
 514      &        uoce=uoce(:,j),voce=voce(:,j),                           & !input
 515      &        qke=QKE(:,:,j),qke_adv=qke_adv(:,:,j),                   & !output
 516      &        sh3d=Sh3d(:,:,j),sm3d=Sm3d(:,:,j),                       & !output
 517      &        nchem=nchem,kdvel=kdvel,ndvel=ndvel,                     & !chem/smoke
 518      &        Chem3d=chem,Vdep=vd,                                     &
 519      &        FRP=frp_mean,EMIS_ANT_NO=emis_ant_no,                    &
 520      &        mix_chem=mix_chem,enh_mix=enh_mix,                       &
 521      &        rrfs_sd=rrfs_sd,smoke_dbg=smoke_dbg,                     & !end chem/smoke
 522      &        tsq=tsq(:,:,j),qsq=qsq(:,:,j),cov=cov(:,:,j),            & !output
 523      &        RUBLTEN=RUBLTEN(:,:,j),RVBLTEN=RVBLTEN(:,:,j),           & !output
 524      &        RTHBLTEN=RTHBLTEN(:,:,j),RQVBLTEN=RQVBLTEN(:,:,j),       & !output
 525      &        RQCBLTEN=rqcblten(:,:,j),RQIBLTEN=rqiblten(:,:,j),       & !output
 526      &        RQNCBLTEN=rqncblten(:,:,j),RQNIBLTEN=rqniblten(:,:,j),   & !output
 527      &        RQNWFABLTEN=RQNWFABLTEN(:,:,j),                          & !output
 528      &        RQNIFABLTEN=RQNIFABLTEN(:,:,j),                          & !output
 529      &        RQNBCABLTEN=RQNBCABLTEN(:,:,j),                          & !output
 530      &        dozone=r03blten(:,:,j),                                  & !output
 531      &        EXCH_H=exch_h(:,:,j),EXCH_M=exch_m(:,:,j),               & !output
 532      &        pblh=pblh(:,j),KPBL=KPBL(:,j),                           & !output
 533      &        el_pbl=el_pbl(:,:,j),                                    & !output
 534      &        dqke=dqke2d,qWT=qWT2d,qSHEAR=qSHEAR2d,                   & !output
 535      &        qBUOY=qBUOY2d,qDISS=qDISS2d,                             & !output
 536      &        qc_bl=qc_bl2d,qi_bl=qi_bl2d,cldfra_bl=cldfra_bl2d,       & !output
 537      &        bl_mynn_tkeadvect=bl_mynn_tkeadvect,                     & !input parameter
 538      &        tke_budget=tke_budget,                                   & !input parameter
 539      &        bl_mynn_cloudpdf=bl_mynn_cloudpdf,                       & !input parameter
 540      &        bl_mynn_mixlength=bl_mynn_mixlength,                     & !input parameter
 541      &        icloud_bl=icloud_bl,                                     & !input parameter
 542      &        closure=bl_mynn_closure,bl_mynn_edmf=bl_mynn_edmf,       & !input parameter
 543      &        bl_mynn_edmf_mom=bl_mynn_edmf_mom,                       & !input parameter
 544      &        bl_mynn_edmf_tke=bl_mynn_edmf_tke,                       & !input parameter
 545      &        bl_mynn_mixscalars=bl_mynn_mixscalars,                   & !input parameter
 546      &        bl_mynn_output=bl_mynn_output,                           & !input parameter
 547      &        bl_mynn_cloudmix=bl_mynn_cloudmix,                       & !input parameter
 548      &        bl_mynn_mixqt=bl_mynn_mixqt,                             & !input parameter
 549      &        edmf_a=edmf_a2d,edmf_w=edmf_w2d,                         & !output
 550      &        edmf_qt=edmf_qt2d,edmf_thl=edmf_thl2d,                   & !output
 551      &        edmf_ent=edmf_ent2d,edmf_qc=edmf_qc2d,                   & !output
 552      &        sub_thl3D=sub_thl2d,sub_sqv3D=sub_sqv2d,                 & !output
 553      &        det_thl3D=det_thl2d,det_sqv3D=det_sqv2d,                 & !output
 554      &        nupdraft=nupdraft(:,j),maxMF=maxMF(:,j),                 & !output
 555      &        ktop_plume=ktop_plume(:,j),                              & !output
 556      &        spp_pbl=spp_pbl,pattern_spp_pbl=pattern_spp_pbl2d,       & !input
 557      &        RTHRATEN=rthraten(:,:,j),                                & !input
 558      &        FLAG_QI=flag_qi,FLAG_QNI=flag_qni,                       & !input
 559      &        FLAG_QC=flag_qc,FLAG_QNC=flag_qnc,                       & !input
 560      &        FLAG_QNWFA=FLAG_QNWFA,FLAG_QNIFA=FLAG_QNIFA,             & !input
 561      &        FLAG_QNBCA=FLAG_QNBCA,                                   & !input
 562      &        IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde,         & !input
 563      &        IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme,         & !input
 564      &        ITS=its,ITE=itf,JTS=jts,JTE=jtf,KTS=kts,KTE=kte)           !input
 565 !print*,"In mynn wrapper, after bl_mynn_driver"
 566
 567      !- Convert spec hum to mixing ratio:
 568       do k=kts,ktf
 569         do i=its,itf
 570            RQVBLTEN(i,k,j) = RQVBLTEN(i,k,j)/(1.0 - sqv(i,k))
 571            RQCBLTEN(i,k,j) = RQCBLTEN(i,k,j)/(1.0 - sqv(i,k))
 572            RQIBLTEN(i,k,j) = RQIBLTEN(i,k,j)/(1.0 - sqv(i,k))
 573         enddo
 574       enddo
 575
 576      !- Collect 3D ouput:
 577       if (icloud_bl > 0) then
 578          do k=kts,ktf
 579          do i=its,itf
 580             qc_bl(i,k,j)     = qc_bl2d(i,k)
 581             qi_bl(i,k,j)     = qi_bl2d(i,k)
 582             cldfra_bl(i,k,j) = cldfra_bl2d(i,k)
 583          enddo
 584          enddo
 585       endif
 586
 587       if (tke_budget .eq. 1) then
 588          do k=kts,ktf
 589          do i=its,itf
 590             dqke(i,k,j)     = dqke2d(i,k)
 591             qwt(i,k,j)      = qwt2d(i,k)
 592             qshear(i,k,j)   = qshear2d(i,k)
 593             qbuoy(i,k,j)    = qbuoy2d(i,k)
 594             qdiss(i,k,j)    = qdiss2d(i,k)
 595          enddo
 596          enddo
 597       endif
 598
 599       if (bl_mynn_output > 0) then
 600          do k=kts,ktf
 601          do i=its,itf
 602             edmf_a(i,k,j)     = edmf_a2d(i,k)
 603             edmf_w(i,k,j)     = edmf_w2d(i,k)
 604             edmf_qt(i,k,j)    = edmf_qt2d(i,k)
 605             edmf_thl(i,k,j)   = edmf_thl2d(i,k)
 606             edmf_ent(i,k,j)   = edmf_ent2d(i,k)
 607             edmf_qc(i,k,j)    = edmf_qc2d(i,k)
 608             sub_thl3d(i,k,j)  = sub_thl2d(i,k)
 609             sub_sqv3d(i,k,j)  = sub_sqv2d(i,k)
 610             det_thl3d(i,k,j)  = det_thl2d(i,k)
 611             det_sqv3d(i,k,j)  = det_sqv2d(i,k)
 612          enddo
 613          enddo
 614       endif
 615
 616       if (debug) then
 617           print*
 618           print*,"===Finished with mynn_bl_driver; output:"
 619           print*,"T:",t3d(its,1,j),t3d(its,2,j),t3d(its,kte,j)
 620           print*,"TH:",th(its,1,j),th(its,2,j),th(its,kte,j)
 621           print*,"rho:",rho(its,1,j),rho(its,2,j),rho(its,kte,j)
 622           print*,"exner:",exner(its,1,j),exner(its,2,j),exner(its,kte,j)
 623           print*,"p:",p(its,1,j),p(its,2,j),p(its,kte,j)
 624           print*,"dz:",dz(its,1,j),dz(its,2,j),dz(its,kte,j)
 625           print*,"u:",u(its,1,j),u(its,2,j),u(its,kte,j)
 626           print*,"v:",v(its,1,j),v(its,2,j),v(its,kte,j)
 627           print*,"sqv:",sqv(its,1),sqv(its,2),sqv(its,kte)
 628           print*,"sqc:",sqc(its,1),sqc(its,2),sqc(its,kte)
 629           print*,"sqi:",sqi(its,1),sqi(its,2),sqi(its,kte)
 630           print*,"rmol:",rmol(its,j)," ust:",ust(its,j)
 631           print*,"dx(its,j)=",dx(its),"initflag=",initflag
 632           print*,"Thetasurf:",ts(its,j)
 633           print*,"HFX:",hfx(its,j)," qfx",qfx(its,j)
 634           print*,"qsfc:",qsfc(its,j)," ps:",ps(its,j)
 635           print*,"wspd:",wspd(its,j)
 636           print*,"znt:",znt(its,j)," delt=",delt
 637           print*,"im=",ite," kte=",kte
 638           print*,"PBLH=",pblh(its,j)," KPBL=",KPBL(its,j)," xland=",xland(its,j)
 639           print*,"ch=",ch(its,j)
 640           print*,"qke:",qke(its,1,j),qke(its,2,j),qke(its,kte,j)
 641           print*,"el_pbl:",el_pbl(its,1,j),el_pbl(its,2,j),el_pbl(its,kte,j)
 642           print*,"Sh3d:",Sh3d(its,1,j),sh3d(its,2,j),sh3d(its,kte,j)
 643           print*,"exch_h:",exch_h(its,1,j),exch_h(its,2,j),exch_h(its,kte,j)
 644           print*,"exch_m:",exch_m(its,1,j),exch_m(its,2,j),exch_m(its,kte,j)
 645           print*,"max cf_bl:",maxval(cldfra_bl(its,:,j))
 646           print*,"max qc_bl:",maxval(qc_bl(its,:,j))
 647           print*,"dtdt:",rthblten(its,1,j),rthblten(its,2,j),rthblten(its,kte,j)
 648           print*,"dudt:",rublten(its,1,j),rublten(its,2,j),rublten(its,kte,j)
 649           print*,"dvdt:",rvblten(its,1,j),rvblten(its,2,j),rvblten(its,kte,j)
 650           print*,"dqdt:",rqvblten(its,1,j),rqvblten(its,2,j),rqvblten(its,kte,j)
 651           print*,"ktop_plume:",ktop_plume(its,j)," maxmf:",maxmf(its,j)
 652           print*,"nup:",nupdraft(its,j)
 653           print*
 654        endif
 655
 656    enddo  !end j-loop
 657
 658    !Deallocate all temporary interface arrays
 659    if (bl_mynn_output > 0) then
 660       deallocate(edmf_a2d)
 661       deallocate(edmf_w2d)
 662       deallocate(edmf_qt2d)
 663       deallocate(edmf_thl2d)
 664       deallocate(edmf_ent2d)
 665       deallocate(edmf_qc2d)
 666       deallocate(sub_thl2d)
 667       deallocate(sub_sqv2d)
 668       deallocate(det_thl2d)
 669       deallocate(det_sqv2d)
 670    endif
 671    if (tke_budget .eq. 1) then
 672       deallocate(dqke2d)
 673       deallocate(qwt2d)
 674       deallocate(qshear2d)
 675       deallocate(qbuoy2d)
 676       deallocate(qdiss2d)
 677    endif
 678    if (icloud_bl > 0) then
 679       deallocate(qc_bl2d)
 680       deallocate(qi_bl2d)
 681       deallocate(cldfra_bl2d)
 682    endif
 683    if (flag_qc)   deallocate(qc2d)
 684    if (flag_qi)   deallocate(qi2d)
 685    if (flag_qnc)  deallocate(qnc2d)
 686    if (flag_qni)  deallocate(qni2d)
 687    if (flag_qnwfa)deallocate(qnwfa2d)
 688    if (flag_qnifa)deallocate(qnifa2d)
 689    if (flag_qnbca)deallocate(qnbca2d)
 690    if (spp_pbl > 0) then
 691       deallocate(pattern_spp_pbl2d)
 692    endif
 693
 694 !print*,"In mynn wrapper, at end"
 695
 696   CONTAINS
 697
 698 ! ==================================================================
 699   SUBROUTINE moisture_check2(kte, delt, dp, exner, &
 700                              qv, qc, qi, th        )
 701   !
 702   ! If qc < qcmin, qi < qimin, or qv < qvmin happens in any layer,
 703   ! force them to be larger than minimum value by (1) condensating
 704   ! water vapor into liquid or ice, and (2) by transporting water vapor
 705   ! from the very lower layer.
 706   !
 707   ! We then update the final state variables and tendencies associated
 708   ! with this correction. If any condensation happens, update theta/temperature too.
 709   ! Note that (qv,qc,qi,th) are the final state variables after
 710   ! applying corresponding input tendencies and corrective tendencies.
 711
 712     implicit none
 713     integer,  intent(in)     :: kte
 714     real, intent(in)     :: delt
 715     real, dimension(kte), intent(in)     :: dp
 716     real, dimension(kte), intent(in)     :: exner
 717     real, dimension(kte), intent(inout)  :: qv, qc, qi, th
 718     integer   k
 719     real ::  dqc2, dqi2, dqv2, sum, aa, dum
 720     real, parameter :: qvmin1= 1e-8,    & !min at k=1
 721                        qvmin = 1e-20,   & !min above k=1
 722                        qcmin = 0.0,     &
 723                        qimin = 0.0
 724
 725     do k = kte, 1, -1  ! From the top to the surface
 726        dqc2 = max(0.0, qcmin-qc(k)) !qc deficit (>=0)
 727        dqi2 = max(0.0, qimin-qi(k)) !qi deficit (>=0)
 728
 729        !update species
 730        qc(k)  = qc(k)  +  dqc2
 731        qi(k)  = qi(k)  +  dqi2
 732        qv(k)  = qv(k)  -  dqc2 - dqi2
 733        !for theta
 734        !th(k)  = th(k)  +  xlvcp/exner(k)*dqc2 + &
 735        !                   xlscp/exner(k)*dqi2
 736        !for temperature
 737        th(k)  = th(k)  +  xlvcp*dqc2 + &
 738                           xlscp*dqi2
 739
 740        !then fix qv if lending qv made it negative
 741        if (k .eq. 1) then
 742           dqv2   = max(0.0, qvmin1-qv(k)) !qv deficit (>=0)
 743           qv(k)  = qv(k)  + dqv2
 744           qv(k)  = max(qv(k),qvmin1)
 745           dqv2   = 0.0
 746        else
 747           dqv2   = max(0.0, qvmin-qv(k)) !qv deficit (>=0)
 748           qv(k)  = qv(k)  + dqv2
 749           qv(k-1)= qv(k-1)  - dqv2*dp(k)/dp(k-1)
 750           qv(k)  = max(qv(k),qvmin)
 751        endif
 752        qc(k) = max(qc(k),qcmin)
 753        qi(k) = max(qi(k),qimin)
 754     end do
 755
 756     ! Extra moisture used to satisfy 'qv(1)>=qvmin' is proportionally
 757     ! extracted from all the layers that has 'qv > 2*qvmin'. This fully
 758     ! preserves column moisture.
 759     if( dqv2 .gt. 1.e-20 ) then
 760         sum = 0.0
 761         do k = 1, kte
 762            if( qv(k) .gt. 2.0*qvmin ) sum = sum + qv(k)*dp(k)
 763         enddo
 764         aa = dqv2*dp(1)/max(1.e-20,sum)
 765         if( aa .lt. 0.5 ) then
 766             do k = 1, kte
 767                if( qv(k) .gt. 2.0*qvmin ) then
 768                    dum    = aa*qv(k)
 769                    qv(k)  = qv(k) - dum
 770                endif
 771             enddo
 772         else
 773         ! For testing purposes only (not yet found in any output):
 774         !    write(*,*) 'Full moisture conservation is impossible'
 775         endif
 776     endif
 777
 778     return
 779
 780   END SUBROUTINE moisture_check2
 781
 782   END SUBROUTINE mynnedmf_wrapper_run
 783
 784 !###=================================================================
 785
 786 END MODULE module_bl_mynn_wrapper