ufs-dev PR#371

.gitmodules

@@ -9,4 +9,4 @@
 [submodule "physics/SFC_Layer/MYNN/MYNN"]
   path = physics/SFC_Layer/MYNN/MYNN
   url = https://github.com/NCAR/MYNN-SFC
-  branch = main
+  branch = ccpp/dev

physics/CONV/C3/cu_c3_deep.F90

@@ -744,20 +744,18 @@ subroutine cu_c3_deep_run(        &
         do k=kts,ktf
           if(zo_cup(i,k).gt.zkbmax+z1(i))then
             kbmax(i)=k
-            go to 25
+            exit
           endif
         enddo
- 25     continue
 !
 !> - Compute the level where detrainment for downdraft starts (\p kdet)
 !
         do k=kts,ktf
           if(zo_cup(i,k).gt.z_detr+z1(i))then
             kdet(i)=k
-            go to 26
+            exit
           endif
         enddo
- 26     continue
 !
         endif
       enddo
@@ -770,7 +768,7 @@ subroutine cu_c3_deep_run(        &
 !
       start_k22=2
 !$acc parallel loop
-       do 36 i=its,itf
+       do i=its,itf
          if(ierr(i).eq.0)then
             k22(i)=maxloc(heo_cup(i,start_k22:kbmax(i)+2),1)+start_k22-1
             if(k22(i).ge.kbmax(i))then
@@ -783,7 +781,7 @@ subroutine cu_c3_deep_run(        &
              kbcon(i)=0
            endif
          endif
- 36   continue
+       end do
 !$acc end parallel
 
 !
@@ -1013,7 +1011,7 @@ subroutine cu_c3_deep_run(        &
 !$acc end parallel
 
 !$acc kernels
-      do 37 i=its,itf
+      do i=its,itf
          kzdown(i)=0
          if(ierr(i).eq.0)then
             zktop=(zo_cup(i,ktop(i))-z1(i))*.6
@@ -1024,11 +1022,11 @@ subroutine cu_c3_deep_run(        &
               if(zo_cup(i,k).gt.zktop)then
                  kzdown(i)=k
                  kzdown(i)=min(kzdown(i),kstabi(i)-1)  !
-                 go to 37
+                 exit
               endif
               enddo
          endif
- 37   continue
+      end do
 !$acc end kernels
 
 !
@@ -2670,8 +2668,8 @@ subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav,            &
         edtc(i,1)=0.
        enddo
        do kk = kts,ktf-1
-         do 62 i=its,itf
-          if(ierr(i).ne.0)go to 62
+         do i=its,itf
+          if(ierr(i).ne.0) cycle
           if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
              vws(i) = vws(i)+                                        &
               (abs((us(i,kk+1)-us(i,kk))/(z(i,kk+1)-z(i,kk)))        &
@@ -2680,7 +2678,7 @@ subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav,            &
             sdp(i) = sdp(i) + p(i,kk) - p(i,kk+1)
           endif
           if (kk .eq. ktf-1)vshear(i) = 1.e3 * vws(i) / sdp(i)
-   62   continue
+         end do
        end do
       do i=its,itf
          if(ierr(i).eq.0)then
@@ -3632,13 +3630,13 @@ subroutine cup_kbcon(ierrc,cap_inc,iloop_in,k22,kbcon,he_cup,hes_cup, &
 !$acc end kernels
 
 !$acc parallel loop
-       do 27 i=its,itf
+      do i=its,itf
       kbcon(i)=1
 !
 ! reset iloop for mid level convection
       if(cap_max(i).gt.200 .and. imid.eq.1)iloop(i)=5
 !
-      if(ierr(i).ne.0)go to 27
+      if(ierr(i).ne.0) cycle
       start_level(i)=k22(i)
       kbcon(i)=k22(i)+1
       if(iloop(i).eq.5)kbcon(i)=k22(i)
@@ -3654,70 +3652,73 @@ subroutine cup_kbcon(ierrc,cap_inc,iloop_in,k22,kbcon,he_cup,hes_cup, &
         enddo
        !==
 
-      go to 32
- 31   continue
-      kbcon(i)=kbcon(i)+1
-      if(kbcon(i).gt.kbmax(i)+2)then
-         if(iloop(i).ne.4)then
-                ierr(i)=3
+       find_kbcon_loop: do
+
+        hetest=hcot(i,kbcon(i)) !hkb(i) ! he_cup(i,k22(i))
+        if(hetest.lt.hes_cup(i,kbcon(i)))then
+          kbcon(i)=kbcon(i)+1
+          if(kbcon(i).gt.kbmax(i)+2)then
+             if(iloop(i).ne.4)then
+                    ierr(i)=3
 #ifndef _OPENACC
-                ierrc(i)="could not find reasonable kbcon in cup_kbcon"
+                    ierrc(i)="could not find reasonable kbcon in cup_kbcon"
 #endif
-         endif
-        go to 27
-      endif
- 32   continue
-      hetest=hcot(i,kbcon(i)) !hkb(i) ! he_cup(i,k22(i))
-      if(hetest.lt.hes_cup(i,kbcon(i)))then
-        go to 31
-      endif
+             endif
+             exit find_kbcon_loop
+          endif
+          cycle find_kbcon_loop
+        endif
 
-!     cloud base pressure and max moist static energy pressure
-!     i.e., the depth (in mb) of the layer of negative buoyancy
-      if(kbcon(i)-k22(i).eq.1)go to 27
-      if(iloop(i).eq.5 .and. (kbcon(i)-k22(i)).le.2)go to 27
-      pbcdif=-p_cup(i,kbcon(i))+p_cup(i,k22(i))
-      plus=max(25.,cap_max(i)-float(iloop(i)-1)*cap_inc(i))
-      if(iloop(i).eq.4)plus=cap_max(i)
+!      cloud base pressure and max moist static energy pressure
+!      i.e., the depth (in mb) of the layer of negative buoyancy
+        if(kbcon(i)-k22(i).eq.1) exit find_kbcon_loop
+        if(iloop(i).eq.5 .and. (kbcon(i)-k22(i)).le.2) exit find_kbcon_loop
+
+        pbcdif=-p_cup(i,kbcon(i))+p_cup(i,k22(i))
+        plus=max(25.,cap_max(i)-float(iloop(i)-1)*cap_inc(i))
+        if(iloop(i).eq.4)plus=cap_max(i)
 !
 ! for shallow convection, if cap_max is greater than 25, it is the pressure at pbltop
-      if(iloop(i).eq.5)plus=150.
-        if(iloop(i).eq.5.and.cap_max(i).gt.200)pbcdif=-p_cup(i,kbcon(i))+cap_max(i)
-      if(pbcdif.le.plus)then
-        go to 27
-      elseif(pbcdif.gt.plus)then
-        k22(i)=k22(i)+1
-        kbcon(i)=k22(i)+1
-!==     since k22 has be changed, hkb has to be re-calculated
-        x_add = xlv*zqexec(i)+cp*ztexec(i)
-        call get_cloud_bc(kte,he_cup (i,1:kte),hkb (i),k22(i),x_add)
-
-        start_level(i)=k22(i)
-!        if(iloop_in.eq.5)start_level(i)=kbcon(i)
-        hcot(i,1:start_level(i)) = hkb(i)
+        if(iloop(i).eq.5)plus=150.
+          if(iloop(i).eq.5.and.cap_max(i).gt.200)pbcdif=-p_cup(i,kbcon(i))+cap_max(i)
+
+        if(pbcdif.le.plus)then
+          exit find_kbcon_loop
+        elseif(pbcdif.gt.plus)then
+          k22(i)=k22(i)+1
+          kbcon(i)=k22(i)+1
+!==      since k22 has be changed, hkb has to be re-calculated
+          x_add = xlv*zqexec(i)+cp*ztexec(i)
+          call get_cloud_bc(kte,he_cup (i,1:kte),hkb (i),k22(i),x_add)
+
+          start_level(i)=k22(i)
+!          if(iloop_in.eq.5)start_level(i)=kbcon(i)
+          hcot(i,1:start_level(i)) = hkb(i)
 !$acc loop seq
-        do k=start_level(i)+1,kbmax(i)+3
-           dz=z_cup(i,k)-z_cup(i,k-1)
+          do k=start_level(i)+1,kbmax(i)+3
+             dz=z_cup(i,k)-z_cup(i,k-1)
 
-           hcot(i,k)= ( (1.-0.5*entr_rate(i)*dz)*hcot(i,k-1)   &
-                         + entr_rate(i)*dz*heo(i,k-1)       )/ &
-                      (1.+0.5*entr_rate(i)*dz)
-        enddo
-       !==
+             hcot(i,k)= ( (1.-0.5*entr_rate(i)*dz)*hcot(i,k-1)   &
+                           + entr_rate(i)*dz*heo(i,k-1)        )/ &
+                        (1.+0.5*entr_rate(i)*dz)
+          enddo
+          !==
 
-        if(iloop(i).eq.5)kbcon(i)=k22(i)
-        if(kbcon(i).gt.kbmax(i)+2)then
-            if(iloop(i).ne.4)then
-                ierr(i)=3
+          if(iloop(i).eq.5)kbcon(i)=k22(i)
+          if(kbcon(i).gt.kbmax(i)+2)then
+              if(iloop(i).ne.4)then
+                  ierr(i)=3
 #ifndef _OPENACC
-                ierrc(i)="could not find reasonable kbcon in cup_kbcon"
+                  ierrc(i)="could not find reasonable kbcon in cup_kbcon"
 #endif
-            endif
-            go to 27
+              endif
+              exit find_kbcon_loop
+          endif
+          cycle find_kbcon_loop
         endif
-        go to 32
-      endif
- 27   continue
+      end do find_kbcon_loop
+
+      end do
  !$acc end parallel
 
    end subroutine cup_kbcon
@@ -3946,7 +3947,7 @@ subroutine neg_check(name,j,dt,q,outq,outt,outu,outv,                      &
         ,intent (in  )                   ::                                   &
         dt
      real(kind=kind_phys) :: names,scalef,thresh,qmem,qmemf,qmem2,qtest,qmem1
-     integer :: icheck
+     integer :: icheck, i, k
 !
 ! first do check on vertical heating rate
 !
@@ -4815,17 +4816,18 @@ subroutine rates_up_pdf(rand_vmas,ipr,name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo
         enddo
         ktopdby(i)=maxloc(dby(:),1)
         kklev=maxloc(dbm(:),1)
+
+        kfinalzu=ktf-2
+        ktop(i)=kfinalzu
 !$acc loop seq
         do k=maxloc(dby(:),1)+1,ktf-2
           if(dby(k).lt.dbythresh*maxval(dby))then
               kfinalzu=k  - 1
               ktop(i)=kfinalzu
-              go to 412
+              exit
           endif
         enddo
-        kfinalzu=ktf-2
-        ktop(i)=kfinalzu
-412     continue
+
         ktop(i)=ktopdby(i) ! HCB
         kklev=min(kklev+3,ktop(i)-2)
 !
@@ -5774,10 +5776,9 @@ subroutine get_cloud_top(name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo,heso_cup,z_c
               kfinalzu = k - 1
               ktop(i)  = kfinalzu
               !print*,'hco4=',k,kfinalzu,ktop(i),kbcon(i)+1;call flush(6)
-              go to 412
+              exit
           endif
         enddo
-        412    continue
        else
          do k=start_level(i)+1,ktf-2
           !~ print*,'hco31=',k,dby(k),dbythresh*maxval(dby)

physics/CONV/C3/cu_c3_driver.F90

@@ -80,7 +80,7 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                fhour,fh_dfi_radar,ix_dfi_radar,num_dfi_radar,cap_suppress,      &
                dfi_radar_max_intervals,ldiag3d,qci_conv,do_cap_suppress,        &
                sigmaout,maxupmf,maxMF,do_mynnedmf,ichoice_in,ichoicem_in,       &
-               ichoice_s_in,errmsg,errflg)
+               ichoice_s_in,ten_t,ten_u,ten_v,ten_q,dcliw,dclcw,errmsg,errflg)
 !-------------------------------------------------------------
       implicit none
       integer, parameter :: maxiens=1
@@ -120,11 +120,11 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
    real(kind=kind_phys),  dimension( : , : ), intent(in    ), optional :: forcet,forceqv_spechum
    real(kind=kind_phys),  dimension( : , : ), intent(in    ) :: w,phil,delp
    real(kind=kind_phys), dimension ( : , : ), intent(in    ), optional :: sigmain,qmicro
-   real(kind=kind_phys),  dimension( : , : ), intent(inout ) :: t,us,vs
+   real(kind=kind_phys),  dimension( : , : ), intent(in    ) :: t,us,vs
    real(kind=kind_phys),  dimension( : , : ), intent(inout ), optional :: qci_conv
    real(kind=kind_phys),  dimension( : , : ), intent(out   ) :: cnvw_moist,cnvc
    real(kind=kind_phys), dimension ( : , : ), intent(out   ), optional :: sigmaout
-   real(kind=kind_phys),  dimension( : , : ), intent(inout ) :: cliw, clcw
+   real(kind=kind_phys),  dimension( : , : ), intent(in    ) :: cliw, clcw
    real(kind=kind_phys), dimension ( : , : , :), intent(in    ) :: tmf
 !$acc declare copyin(forcet,forceqv_spechum,w,phil)
 !$acc declare copy(t,us,vs,qci_conv,cliw, clcw)
@@ -158,7 +158,7 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
 !$acc declare copyout(ud_mf,dd_mf,dt_mf,raincv,cld1d)
    ! Specific humidity from FV3
    real(kind=kind_phys), dimension (:,:), intent(in) :: qv2di_spechum
-   real(kind=kind_phys), dimension (:,:), intent(inout) :: qv_spechum
+   real(kind=kind_phys), dimension (:,:), intent(in) :: qv_spechum
    real(kind=kind_phys), dimension (:), intent(inout), optional :: aod_gf
 !$acc declare copyin(qv2di_spechum) copy(qv_spechum,aod_gf)
    ! Local water vapor mixing ratios and cloud water mixing ratios
@@ -172,7 +172,9 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
    integer, intent(in   ) :: imfshalcnv
    integer, dimension(:), intent(inout), optional :: cactiv,cactiv_m
 !$acc declare copy(cactiv,cactiv_m)
-
+   real(kind_phys), dimension(:,:), intent(out) :: ten_t, ten_u, ten_v, dcliw, dclcw
+   real(kind_phys), dimension(:,:,:), intent(out) :: ten_q
+
    character(len=*), intent(out) :: errmsg
    integer,          intent(out) :: errflg
 
@@ -225,7 +227,7 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
 ! omega (omeg), windspeed (us,vs), and a flag (ierr) to turn off
 ! convection for this call only and at that particular gridpoint
 !
-   real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi,clw_ten
+   real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi,clw_ten,new_qv_spechum,new_cliw,new_clcw
    real(kind=kind_phys), dimension (im,km) :: tn,qo,tshall,qshall,dz8w,omeg
    real(kind=kind_phys), dimension (im)    :: z1,psur,cuten,cutens,cutenm
    real(kind=kind_phys), dimension (im)    :: umean,vmean,pmean,mc_thresh
@@ -268,7 +270,16 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
   ! initialize ccpp error handling variables
      errmsg = ''
      errflg = 0
-
+
+     ten_t = 0.0
+     ten_u = 0.0
+     ten_v = 0.0
+     ten_q = 0.0
+     dcliw = 0.0
+     dclcw = 0.0
+     new_clcw = clcw
+     new_cliw = cliw
+
      ichoice   = ichoice_in
      ichoicem  = ichoicem_in
      ichoice_s = ichoice_s_in
@@ -956,11 +967,12 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                cnvw(i,k)=cnvwt(i,k)*xmb(i)*dt+cnvwts(i,k)*xmbs(i)*dt+cnvwtm(i,k)*xmbm(i)*dt
                ud_mf(i,k)=cuten(i)*zu(i,k)*xmb(i)*dt
                dd_mf(i,k)=cuten(i)*zd(i,k)*edt(i)*xmb(i)*dt
-               t(i,k)=t(i,k)+dt*(cutens(i)*outts(i,k)+cutenm(i)*outtm(i,k)+outt(i,k)*cuten(i))
+
+               ten_t(i,k) = cutens(i)*outts(i,k)+cutenm(i)*outtm(i,k)+outt(i,k)*cuten(i)
                qv(i,k)=max(1.e-16,qv(i,k)+dt*(cutens(i)*outqs(i,k)+cutenm(i)*outqm(i,k)+outq(i,k)*cuten(i)))
                gdc(i,k,7)=sqrt(us(i,k)**2 +vs(i,k)**2)
-               us(i,k)=us(i,k)+outu(i,k)*cuten(i)*dt +outum(i,k)*cutenm(i)*dt +outus(i,k)*cutens(i)*dt
-               vs(i,k)=vs(i,k)+outv(i,k)*cuten(i)*dt +outvm(i,k)*cutenm(i)*dt +outvs(i,k)*cutens(i)*dt
+               ten_u(i,k) = outu(i,k)*cuten(i) +outum(i,k)*cutenm(i) +outus(i,k)*cutens(i)
+               ten_v(i,k) = outv(i,k)*cuten(i) +outvm(i,k)*cutenm(i) +outvs(i,k)*cutens(i)
 
                gdc(i,k,1)= max(0.,tun_rad_shall(i)*cupclws(i,k)*cutens(i))      ! my mod
                !gdc2(i,k,1)=max(0.,tun_rad_deep(i)*(cupclwm(i,k)*cutenm(i)+cupclw(i,k)*cuten(i)))
@@ -1005,10 +1017,13 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                          )
                tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf))
                if (clcw(i,k) .gt. -999.0) then
-                cliw(i,k) = max(0.,cliw(i,k) + tem * tem1)            ! ice
-                clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1))       ! water
+                new_cliw(i,k) = max(0.,cliw(i,k) + tem * tem1)            ! ice
+                new_clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1))       ! water
+                dcliw(i,k) = (new_cliw(i,k) - cliw(i,k))/dt
+                dclcw(i,k) = (new_clcw(i,k) - clcw(i,k))/dt
                else
-                cliw(i,k) = max(0.,cliw(i,k) + tem)
+                new_cliw(i,k) = max(0.,cliw(i,k) + tem)
+                dcliw(i,k) = (new_cliw(i,k) - cliw(i,k))/dt
                endif
 
              enddo
@@ -1076,9 +1091,14 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
 ! Scale dry mixing ratios for water wapor and cloud water to specific humidy / moist mixing ratios
 !
 !$acc kernels
-        qv_spechum = qv/(1.0_kind_phys+qv)
+        new_qv_spechum = qv/(1.0_kind_phys+qv)
         cnvw_moist = cnvw/(1.0_kind_phys+qv)
 !$acc end kernels
+        do i=its,ite
+          do k=kts,kte
+            ten_q(i,k,ntqv) = (new_qv_spechum(i,k) - qv_spechum(i,k))/dt
+          end do
+        end do
 !
 ! Diagnostic tendency updates
 !