LCOV - code coverage report
Current view: | top level - atmos_phys/schemes/dry_adiabatic_adjust - dadadj.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 50 | 89 | 56.2 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 2 | 2 | 100.0 %

              Line data    Source code
           1             : module dadadj
           2             :   !=======================================================================
           3             :   ! GFDL style dry adiabatic adjustment
           4             :   !
           5             :   ! Method:
           6             :   ! if stratification is unstable, adjustment to the dry adiabatic lapse
           7             :   ! rate is forced subject to the condition that enthalpy is conserved.
           8             :   !=======================================================================
           9             :
          10             :   use ccpp_kinds, only:  kind_phys
          11             :
          12             :   implicit none
          13             :   private
          14             :   save
          15             :
          16             :   public :: dadadj_init ! init routine
          17             :   public :: dadadj_run  ! main routine
          18             :
          19             :   integer  :: nlvdry  ! number of layers from top of model to apply the adjustment
          20             :   integer  :: niter   ! number of iterations for convergence
          21             :
          22             : CONTAINS
          23             :
          24             :   !> \section arg_table_dadadj_init Argument Table
          25             :   !! \htmlinclude dadadj_init.html
          26        1536 :   subroutine dadadj_init(dadadj_nlvdry, dadadj_niter, nz, errmsg, errflg)
          27             :     !------------------------------------------------
          28             :     !   Input / output parameters
          29             :     !------------------------------------------------
          30             :     integer,            intent(in)  :: dadadj_nlvdry
          31             :     integer,            intent(in)  :: dadadj_niter
          32             :     integer,            intent(in)  :: nz
          33             :     character(len=512), intent(out) :: errmsg
          34             :     integer,            intent(out) :: errflg
          35             :
          36        1536 :     errmsg = ''
          37        1536 :     errflg = 0
          38             :
          39        1536 :     if (dadadj_nlvdry >= nz .or. dadadj_nlvdry < 0) then
          40           0 :        errflg = 1
          41           0 :        write(errmsg,*) 'dadadj_init: dadadj_nlvdry=',dadadj_nlvdry,' but must be less than the number of vertical levels ',&
          42           0 :        '(',nz,'), and must be a positive integer.'
          43             :     end if
          44             :
          45        1536 :     nlvdry = dadadj_nlvdry
          46        1536 :     niter = dadadj_niter
          47             :
          48        1536 :   end subroutine dadadj_init
          49             :
          50             :   !> \section arg_table_dadadj_run Argument Table
          51             :   !! \htmlinclude dadadj_run.html
          52     1495368 :   subroutine dadadj_run( &
          53     2990736 :        ncol, nz, dt, pmid, pint, pdel, state_t, state_q, cappa, cpair, s_tend, &
          54     4486104 :        q_tend, dadpdf, scheme_name, errmsg, errflg)
          55             :
          56             :     !------------------------------------------------
          57             :     !   Input / output parameters
          58             :     !------------------------------------------------
          59             :     integer,  intent(in)        :: ncol  ! number of atmospheric columns
          60             :     integer,  intent(in)        :: nz    ! number of atmospheric levels
          61             :     real(kind_phys), intent(in) :: dt          ! physics timestep
          62             :     real(kind_phys), intent(in) :: pmid(:,:)   ! pressure at model levels
          63             :     real(kind_phys), intent(in) :: pint(:,:)   ! pressure at model interfaces
          64             :     real(kind_phys), intent(in) :: pdel(:,:)   ! vertical delta-p
          65             :     real(kind_phys), intent(in) :: cappa(:,:)  ! variable Kappa
          66             :     real(kind_phys), intent(in) :: cpair(:,:)  ! heat capacity of air
          67             :     real(kind_phys), intent(in) :: state_t(:,:)   ! temperature (K)
          68             :     real(kind_phys), intent(in) :: state_q(:,:)   ! specific humidity
          69             :     real(kind_phys), intent(out) :: s_tend(:,:)   ! dry air enthalpy tendency
          70             :     real(kind_phys), intent(out) :: q_tend(:,:)   ! specific humidity tendency
          71             :     real(kind_phys), intent(out) :: dadpdf(:,:)  ! PDF of where adjustments happened
          72             :
          73             :     character(len=64),  intent(out) :: scheme_name
          74             :     character(len=512), intent(out) :: errmsg
          75             :     integer,            intent(out) :: errflg
          76             :
          77             :     !------------------------------------------------
          78             :     !   Local variables
          79             :     !------------------------------------------------
          80             :
          81             :     integer                      :: i,k      ! longitude, level indices
          82             :     integer                      :: jiter    ! iteration index
          83     1495368 :     real(kind_phys), allocatable :: c1dad(:) ! intermediate constant
          84     1495368 :     real(kind_phys), allocatable :: c2dad(:) ! intermediate constant
          85     1495368 :     real(kind_phys), allocatable :: c3dad(:) ! intermediate constant
          86     1495368 :     real(kind_phys), allocatable :: c4dad(:) ! intermediate constant
          87             :     real(kind_phys) :: gammad    ! dry adiabatic lapse rate (deg/Pa)
          88             :     real(kind_phys) :: zeps      ! convergence criterion (deg/Pa)
          89             :     real(kind_phys) :: rdenom    ! reciprocal of denominator of expression
          90             :     real(kind_phys) :: dtdp      ! delta-t/delta-p
          91             :     real(kind_phys) :: zepsdp    ! zeps*delta-p
          92             :     real(kind_phys) :: zgamma    ! intermediate constant
          93             :     real(kind_phys) :: qave      ! mean q between levels
          94             :     real(kind_phys) :: cappaint  ! Kappa at level intefaces
          95     2990736 :     real(kind_phys) :: t(ncol,nz)
          96     2990736 :     real(kind_phys) :: q(ncol,nz)
          97             :
          98             :     logical :: ilconv          ! .TRUE. ==> convergence was attained
          99     1495368 :     logical :: dodad(ncol)     ! .TRUE. ==> do dry adjustment
         100             :
         101             :     !-----------------------------------------------------------------------
         102             :
         103     1495368 :     zeps = 2.0e-5_kind_phys           ! set convergence criteria
         104     1495368 :     errmsg = ''
         105     1495368 :     errflg = 0
         106     1495368 :     scheme_name = 'DADADJ'
         107             :
         108     4486104 :     allocate(c1dad(nlvdry), stat=errflg)
         109     1495368 :     if (errflg /= 0) then
         110           0 :        errmsg = trim(scheme_name)//': Allocate of c1dad(nlvdry) failed'
         111           0 :        return
         112             :     end if
         113     4486104 :     allocate(c2dad(nlvdry), stat=errflg)
         114     1495368 :     if (errflg /= 0) then
         115           0 :        errmsg = trim(scheme_name)//': Allocate of c2dad(nlvdry) failed'
         116           0 :        return
         117             :     end if
         118     4486104 :     allocate(c3dad(nlvdry), stat=errflg)
         119     1495368 :     if (errflg /= 0) then
         120           0 :        errmsg = trim(scheme_name)//': Allocate of c3dad(nlvdry) failed'
         121           0 :        return
         122             :     end if
         123     4486104 :     allocate(c4dad(nlvdry), stat=errflg)
         124     1495368 :     if (errflg /= 0) then
         125           0 :        errmsg = trim(scheme_name)//': Allocate of c4dad(nlvdry) failed'
         126           0 :        return
         127             :     end if
         128             :
         129   650693736 :     t = state_t
         130   650693736 :     q = state_q
         131             :
         132             :     ! Find gridpoints with unstable stratification
         133             :
         134    24969168 :     do i = 1, ncol
         135    23473800 :        cappaint = 0.5_kind_phys*(cappa(i,2) + cappa(i,1))
         136    23473800 :        gammad = cappaint*0.5_kind_phys*(t(i,2) + t(i,1))/pint(i,2)
         137    23473800 :        dtdp = (t(i,2) - t(i,1))/(pmid(i,2) - pmid(i,1))
         138    24969168 :        dodad(i) = (dtdp + zeps) > gammad
         139             :     end do
         140             :
         141   650693736 :     dadpdf(:ncol,:) = 0._kind_phys
         142     4486104 :     do k= 2, nlvdry
         143    51433704 :        do i = 1, ncol
         144    46947600 :          cappaint = 0.5_kind_phys*(cappa(i,k+1) + cappa(i,k))
         145    46947600 :          gammad = cappaint*0.5_kind_phys*(t(i,k+1) + t(i,k))/pint(i,k+1)
         146    46947600 :          dtdp = (t(i,k+1) - t(i,k))/(pmid(i,k+1) - pmid(i,k))
         147    46947600 :          dodad(i) = dodad(i) .or. (dtdp + zeps) > gammad
         148    49938336 :          if ((dtdp + zeps) > gammad) then
         149           0 :             dadpdf(i,k) = 1._kind_phys
         150             :          end if
         151             :       end do
         152             :    end do
         153             :
         154             :    ! Make a dry adiabatic adjustment
         155             :    ! Note: nlvdry ****MUST**** be < pver
         156             :
         157    24969168 :    COL: do i = 1, ncol
         158             :
         159    24969168 :       if (dodad(i)) then
         160             :
         161           0 :          zeps = 2.0e-5_kind_phys
         162             :
         163           0 :          do k = 1, nlvdry
         164           0 :             cappaint = 0.5_kind_phys*(cappa(i,k+1) + cappa(i,k))
         165           0 :             c1dad(k) = cappaint*0.5_kind_phys*(pmid(i,k+1)-pmid(i,k))/pint(i,k+1)
         166           0 :             c2dad(k) = (1._kind_phys - c1dad(k))/(1._kind_phys + c1dad(k))
         167           0 :             rdenom = 1._kind_phys/(pdel(i,k)*c2dad(k) + pdel(i,k+1))
         168           0 :             c3dad(k) = rdenom*pdel(i,k)
         169           0 :             c4dad(k) = rdenom*pdel(i,k+1)
         170             :          end do
         171             :
         172             :          ilconv = .false.
         173             :
         174           0 :          DBLZEP: do while (.not. ilconv)
         175             :
         176           0 :             do jiter = 1, niter
         177             :                ilconv = .true.
         178             :
         179           0 :                do k = 1, nlvdry
         180           0 :                   zepsdp = zeps*(pmid(i,k+1) - pmid(i,k))
         181           0 :                   zgamma = c1dad(k)*(t(i,k) + t(i,k+1))
         182             :
         183           0 :                   if ((t(i,k+1)-t(i,k)) >= (zgamma+zepsdp)) then
         184           0 :                      ilconv = .false.
         185           0 :                      t(i,k+1) = t(i,k)*c3dad(k) + t(i,k+1)*c4dad(k)
         186           0 :                      t(i,k) = c2dad(k)*t(i,k+1)
         187           0 :                      qave = (pdel(i,k+1)*q(i,k+1) + pdel(i,k)*q(i,k))/(pdel(i,k+1)+ pdel(i,k))
         188           0 :                      q(i,k+1) = qave
         189           0 :                      q(i,k) = qave
         190             :                   end if
         191             :
         192             :                end do
         193             :
         194           0 :                if (ilconv) cycle COL ! convergence => next longitude
         195             :
         196             :             end do
         197             :
         198           0 :             zeps = zeps + zeps
         199           0 :             if (zeps > 1.e-4_kind_phys) then
         200           0 :                errflg = i
         201           0 :                write(errmsg,*) trim(scheme_name)//': Convergence failure at column ',i,' zeps > 1.e-4 '// &
         202           0 :                     '(errflg set to failing column index)'
         203           0 :                return                ! error return
         204             :             end if
         205             :          end do DBLZEP
         206             :
         207             :       end if
         208             :
         209             :    end do COL
         210             :
         211   650693736 :    s_tend = (t - state_t)/dt*cpair
         212   650693736 :    q_tend = (q - state_q)/dt
         213             :
         214     1495368 :    deallocate(c1dad, c2dad, c3dad, c4dad)
         215             :
         216     1495368 :  end subroutine dadadj_run
         217             :
         218             : end module dadadj

Generated by: LCOV version 1.14