LCOV - code coverage report
Current view: | top level - utils - hycoef.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 93 | 128 | 72.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 4 | 4 | 100.0 %

              Line data    Source code
           1             : module hycoef
           2             :
           3             : use shr_kind_mod,     only: r8 => shr_kind_r8
           4             : use spmd_utils,       only: masterproc
           5             : use pmgrid,           only: plev, plevp
           6             : use cam_logfile,      only: iulog
           7             : use cam_abortutils,   only: endrun
           8             : use pio,              only: file_desc_t, var_desc_t, &
           9             :                             pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, &
          10             :                             pio_double, pio_def_dim, pio_def_var, &
          11             :                             pio_put_var, pio_get_var, &
          12             :                             pio_seterrorhandling, PIO_BCAST_ERROR, PIO_NOERR
          13             :
          14             : implicit none
          15             : private
          16             : save
          17             :
          18             : !-----------------------------------------------------------------------
          19             : !
          20             : ! Purpose: Hybrid level definitions: p = a*p0 + b*ps
          21             : !          interfaces   p(k) = hyai(k)*ps0 + hybi(k)*ps
          22             : !          midpoints    p(k) = hyam(k)*ps0 + hybm(k)*ps
          23             : !
          24             : ! Note: Module data with a target attribute are targets of pointers in hist_coord_t
          25             : !       objects in the cam_history_support module.  They are associated by the calls
          26             : !       to add_hist_coord and add_vert_coord
          27             : !
          28             : !-----------------------------------------------------------------------
          29             :
          30             : real(r8), public, target :: hyai(plevp) ! ps0 component of hybrid coordinate - interfaces
          31             : real(r8), public, target :: hyam(plev)  ! ps0 component of hybrid coordinate - midpoints
          32             : real(r8), public, target :: hybi(plevp) ! ps component of hybrid coordinate - interfaces
          33             : real(r8), public, target :: hybm(plev)  ! ps component of hybrid coordinate - midpoints
          34             :
          35             : real(r8), public :: etamid(plev)      ! hybrid coordinate - midpoints
          36             :
          37             : real(r8), public :: hybd(plev)        ! difference  in b (hybi) across layers
          38             : real(r8), public :: hypi(plevp)       ! reference pressures at interfaces
          39             : real(r8), public :: hypm(plev)        ! reference pressures at midpoints
          40             : real(r8), public :: hypd(plev)        ! reference pressure layer thickness
          41             : #ifdef planet_mars
          42             : real(r8), public, protected :: ps0 = 6.0e1_r8    ! Base state surface pressure (pascals)
          43             : real(r8), public, protected :: psr = 6.0e1_r8    ! Reference surface pressure (pascals)
          44             : #else
          45             : real(r8), public, protected :: ps0 = 1.0e5_r8    ! Base state surface pressure (pascals)
          46             : real(r8), public, protected :: psr = 1.0e5_r8    ! Reference surface pressure (pascals)
          47             : #endif
          48             : real(r8), target :: alev(plev)    ! level values (hPa) for 'lev' coord
          49             : real(r8), target :: ailev(plevp)  ! interface level values for 'ilev' coord
          50             :
          51             : integer, public :: nprlev       ! number of pure pressure levels at top
          52             :
          53             : public hycoef_init
          54             :
          55             : type(var_desc_t) :: hyam_desc, hyai_desc, hybm_desc, hybi_desc, p0_desc
          56             : public init_restart_hycoef, write_restart_hycoef
          57             :
          58             : !=======================================================================
          59             : contains
          60             : !=======================================================================
          61             :
          62        1536 : subroutine hycoef_init(file, psdry)
          63             :
          64             :    use cam_history_support, only: add_hist_coord, add_vert_coord, formula_terms_t
          65             :
          66             :    !-----------------------------------------------------------------------
          67             :    !
          68             :    ! Purpose:
          69             :    ! Defines the locations of model interfaces from input data in the
          70             :    ! hybrid coordinate scheme.  Actual pressure values of model level
          71             :    ! interfaces are determined elsewhere from the fields set here.
          72             :    !
          73             :    ! Method:
          74             :    ! the following fields are set:
          75             :    ! hyai     fraction of reference pressure used for interface pressures
          76             :    ! hyam     fraction of reference pressure used for midpoint pressures
          77             :    ! hybi     fraction of surface pressure used for interface pressures
          78             :    ! hybm     fraction of surface pressure used for midpoint pressures
          79             :    ! hybd     difference of hybi's
          80             :    ! hypi     reference state interface pressures
          81             :    ! hypm     reference state midpoint pressures
          82             :    ! hypd     reference state layer thicknesses
          83             :    ! hypdln   reference state layer thicknesses (log p)
          84             :    ! hyalph   distance from interface to level (used in integrals)
          85             :    ! prsfac   log pressure extrapolation factor (used to compute psl)
          86             :    !
          87             :    ! Author: B. Boville
          88             :    !
          89             :    !-----------------------------------------------------------------------
          90             :
          91             :    ! arguments
          92             :    type(file_desc_t), intent(inout) :: file
          93             :    logical, optional, intent(in)    :: psdry  ! set true when coordinate is based
          94             :                                            ! on dry surface pressure
          95             :
          96             :    ! local variables
          97             :    integer  :: k        ! Level index
          98             :    logical  :: dry_coord
          99             :    real(r8) :: amean, bmean, atest, btest, eps
         100             :    type(formula_terms_t) :: formula_terms ! For the 'lev' and 'ilev' coords
         101             :    !-----------------------------------------------------------------------
         102             :
         103             :    ! check for dry pressure coordinate (default is moist)
         104        1536 :    dry_coord = .false.
         105        1536 :    if (present(psdry)) dry_coord = psdry
         106             :
         107             :    ! read hybrid coeficients
         108        1536 :    call hycoef_read(file)
         109             :
         110             :    ! Set layer locations
         111        1536 :    nprlev = 0
         112      144384 :    do k=1,plev
         113             :
         114             :       ! Interfaces. Set nprlev to the interface above, the first time a
         115             :       ! nonzero surface pressure contribution is found. "nprlev"
         116             :       ! identifies the lowest pure pressure interface.
         117             :
         118      144384 :       if (nprlev==0 .and. hybi(k).ne.0.0_r8) nprlev = k - 1
         119             :    end do
         120             :
         121             :    ! Set nprlev if no nonzero b's have been found. All interfaces are
         122             :    ! pure pressure. A pure pressure model requires other changes as well.
         123        1536 :    if (nprlev==0) nprlev = plev + 2
         124             :
         125             :    ! Set delta sigma part of layer thickness and reference state midpoint
         126             :    ! pressures
         127      144384 :    do k=1,plev
         128      142848 :       hybd(k) = hybi(k+1) - hybi(k)
         129      142848 :       hypm(k) = hyam(k)*ps0 + hybm(k)*psr
         130      144384 :       etamid(k) = hyam(k) + hybm(k)
         131             :    end do
         132             :
         133             :    ! Reference state interface pressures
         134      145920 :    do k=1,plevp
         135      145920 :       hypi(k) = hyai(k)*ps0 + hybi(k)*psr
         136             :    end do
         137             :
         138             :    ! Reference state layer thicknesses
         139      144384 :    do k=1,plev
         140      144384 :       hypd(k) = hypi(k+1) - hypi(k)
         141             :    end do
         142             :
         143             :    ! Test that A's and B's at full levels are arithmetic means of A's and
         144             :    ! B's at interfaces
         145        1536 :    eps    = 1.e-05_r8
         146      144384 :    do k = 1,plev
         147      142848 :       amean = ( hyai(k+1) + hyai(k) )*0.5_r8
         148      142848 :       bmean = ( hybi(k+1) + hybi(k) )*0.5_r8
         149      142848 :       if(amean == 0._r8 .and. hyam(k) == 0._r8) then
         150        1536 :          atest = 0._r8
         151             :       else
         152      141312 :          atest = abs( amean - hyam(k) )/ ( 0.5_r8*( abs(amean + hyam(k)) ) )
         153             :       endif
         154      142848 :       if(bmean == 0._r8 .and. hybm(k) == 0._r8) then
         155       66048 :          btest = 0._r8
         156             :       else
         157       76800 :          btest = abs( bmean - hybm(k) )/ ( 0.5_r8*( abs(bmean + hybm(k)) ) )
         158             :       endif
         159      142848 :       if (atest > eps) then
         160           0 :          if (masterproc) then
         161           0 :             write(iulog,9850)
         162           0 :             write(iulog,*)'k,atest,eps=',k,atest,eps
         163             :          end if
         164             :       end if
         165             :
         166      144384 :       if (btest > eps) then
         167           0 :          if (masterproc) then
         168           0 :             write(iulog,9850)
         169           0 :             write(iulog,*)'k,btest,eps=',k,btest,eps
         170             :          end if
         171             :       end if
         172             :    end do
         173             :
         174             :    ! Add the information for the 'lev' and 'ilev' mdim history coordinates
         175             :    !
         176             :    ! The hybrid coordinate used by the SE dycore is based on a dry surface
         177             :    ! pressure.  Hence it is the dry pressure rather than actual pressure
         178             :    ! that is computed by the formula_terms attribute.  This coordinate is
         179             :    ! not described by the formula
         180             :    ! atmosphere_hybrid_sigma_pressure_coordinate since the formula
         181             :    ! associated with that name uses actual pressure values.  Furthermore,
         182             :    ! the actual pressure field cannot be reconstructed from the hybrid
         183             :    ! coefficients and the surface pressure field.  Hence in the case of a
         184             :    ! dry coordinate we add neither the standard_name nor the formula_terms
         185             :    ! attributes to the lev and ilev coordinates.
         186             :
         187             :    ! 0.01 converts Pascals to millibars
         188      144384 :    alev(:plev) = 0.01_r8*ps0*(hyam(:plev) + hybm(:plev))
         189      145920 :    ailev(:plevp) = 0.01_r8*ps0*(hyai(:plevp) + hybi(:plevp))
         190             :
         191        1536 :    if (dry_coord) then
         192             :       call add_vert_coord('lev', plev,                                       &
         193             :          'hybrid level at midpoints (1000*(A+B))', 'hPa', alev,              &
         194        1536 :          positive='down')
         195             :       call add_hist_coord('hyam', plev, &
         196        1536 :          'hybrid A coefficient at layer midpoints', '1', hyam, dimname='lev')
         197             :       call add_hist_coord('hybm', plev, &
         198        1536 :          'hybrid B coefficient at layer midpoints', '1', hybm, dimname='lev')
         199             :    else
         200             :
         201           0 :       formula_terms%a_name       =  'hyam'
         202           0 :       formula_terms%a_long_name  =  'hybrid A coefficient at layer midpoints'
         203           0 :       formula_terms%a_values     => hyam
         204           0 :       formula_terms%b_name       =  'hybm'
         205           0 :       formula_terms%b_long_name  =  'hybrid B coefficient at layer midpoints'
         206           0 :       formula_terms%b_values     => hybm
         207           0 :       formula_terms%p0_name      =  'P0'
         208           0 :       formula_terms%p0_long_name = 'reference pressure'
         209           0 :       formula_terms%p0_units     =  'Pa'
         210           0 :       formula_terms%p0_value     =  ps0
         211           0 :       formula_terms%ps_name      =  'PS'
         212             :
         213             :       call add_vert_coord('lev', plev,                                       &
         214             :          'hybrid level at midpoints (1000*(A+B))', 'hPa', alev,              &
         215             :          positive='down',                                                    &
         216             :          standard_name='atmosphere_hybrid_sigma_pressure_coordinate',        &
         217           0 :          formula_terms=formula_terms)
         218             :    end if
         219             :
         220        1536 :    if (dry_coord) then
         221             :       call add_vert_coord('ilev', plevp,                                     &
         222             :          'hybrid level at interfaces (1000*(A+B))', 'hPa', ailev,            &
         223        1536 :          positive='down')
         224             :       call add_hist_coord('hyai', plevp, &
         225        1536 :          'hybrid A coefficient at layer interfaces', '1', hyai, dimname='ilev')
         226             :       call add_hist_coord('hybi', plevp, &
         227        1536 :          'hybrid B coefficient at layer interfaces', '1', hybi, dimname='ilev')
         228             :    else
         229           0 :       formula_terms%a_name       =  'hyai'
         230           0 :       formula_terms%a_long_name  =  'hybrid A coefficient at layer interfaces'
         231           0 :       formula_terms%a_values     => hyai
         232           0 :       formula_terms%b_name       =  'hybi'
         233           0 :       formula_terms%b_long_name  =  'hybrid B coefficient at layer interfaces'
         234           0 :       formula_terms%b_values     => hybi
         235           0 :       formula_terms%p0_name      =  'P0'
         236           0 :       formula_terms%p0_long_name = 'reference pressure'
         237           0 :       formula_terms%p0_units     =  'Pa'
         238           0 :       formula_terms%p0_value     =  ps0
         239           0 :       formula_terms%ps_name      =  'PS'
         240             :
         241             :       call add_vert_coord('ilev', plevp,                                     &
         242             :          'hybrid level at interfaces (1000*(A+B))', 'hPa', ailev,            &
         243             :          positive='down',                                                    &
         244             :          standard_name='atmosphere_hybrid_sigma_pressure_coordinate',        &
         245           0 :          formula_terms=formula_terms)
         246             :    end if
         247             :
         248        1536 :    if (masterproc) then
         249           2 :       write(iulog,'(a)')' Layer Locations (*1000) '
         250         188 :       do k=1,plev
         251         186 :          write(iulog,9800)k,hyai(k),hybi(k),hyai(k)+hybi(k)
         252         188 :          write(iulog,9810) hyam(k), hybm(k), hyam(k)+hybm(k)
         253             :       end do
         254             :
         255           2 :       write(iulog,9800)plevp,hyai(plevp),hybi(plevp),hyai(plevp)+hybi(plevp)
         256           2 :       write(iulog,9820)
         257         188 :       do k=1,plev
         258         186 :          write(iulog,9830) k, hypi(k)
         259         188 :          write(iulog,9840) hypm(k), hypd(k)
         260             :       end do
         261           2 :       write(iulog,9830) plevp, hypi(plevp)
         262             :     end if
         263             :
         264             : 9800 format( 1x, i3, 3p, 3(f10.4,10x) )
         265             : 9810 format( 1x, 3x, 3p, 3(10x,f10.4) )
         266             : 9820 format(1x,'reference pressures (Pa)')
         267             : 9830 format(1x,i3,f15.4)
         268             : 9840 format(1x,3x,15x,2f15.4)
         269             : 9850 format('HYCOEF: A and/or B vertical level coefficients at full',/, &
         270             :          ' levels are not the arithmetic mean of half-level values')
         271             :
         272        3072 : end subroutine hycoef_init
         273             :
         274             : !=======================================================================
         275             :
         276        1536 : subroutine init_restart_hycoef(File, vdimids)
         277             :
         278             :    type(file_desc_t), intent(inout) :: File
         279             :    integer,           intent(out)   :: vdimids(:)
         280             :
         281             :    ! PIO traps errors internally, no need to check ierr
         282             :
         283             :    integer :: ierr
         284             :
         285        1536 :    ierr = PIO_Def_Dim(File, 'lev',  plev,  vdimids(1))
         286        1536 :    ierr = PIO_Def_Dim(File, 'ilev', plevp, vdimids(2))
         287             :
         288        1536 :    ierr = pio_def_var(File, 'hyai', pio_double, vdimids(2:2), hyai_desc)
         289        1536 :    ierr = pio_def_var(File, 'hyam', pio_double, vdimids(1:1), hyam_desc)
         290        1536 :    ierr = pio_def_var(File, 'hybi', pio_double, vdimids(2:2), hybi_desc)
         291        1536 :    ierr = pio_def_var(File, 'hybm', pio_double, vdimids(1:1), hybm_desc)
         292             :
         293        1536 :    ierr = pio_def_var(File, 'P0', pio_double, p0_desc)
         294             :
         295        1536 : end subroutine init_restart_hycoef
         296             :
         297             : !=======================================================================
         298             :
         299        1536 : subroutine write_restart_hycoef(file)
         300             :
         301             :    type(file_desc_t), intent(inout) :: File
         302             :
         303             :    ! PIO traps errors internally, no need to check ierr
         304             :
         305             :    integer :: ierr
         306             :
         307        1536 :    ierr = pio_put_var(File, hyai_desc, hyai)
         308        1536 :    ierr = pio_put_var(File, hyam_desc, hyam)
         309        1536 :    ierr = pio_put_var(File, hybi_desc, hybi)
         310        1536 :    ierr = pio_put_var(File, hybm_desc, hybm)
         311             :
         312        1536 :    ierr = pio_put_var(File, p0_desc,   ps0)
         313             :
         314        1536 : end subroutine write_restart_hycoef
         315             :
         316             : !=======================================================================
         317             :
         318        1536 : subroutine hycoef_read(File)
         319             :
         320             :    ! This code is used both for initial and restart reading.
         321             :
         322             :    type(file_desc_t), intent(inout) :: File
         323             :
         324             :    integer :: flev, filev, lev_dimid, ierr
         325             :    integer :: pio_errtype
         326             :
         327             :    type(var_desc_t) :: p0_desc
         328             :
         329             :    character(len=*), parameter :: routine = 'hycoef_read'
         330             :    !----------------------------------------------------------------------------
         331             :
         332             :    ! PIO traps errors internally, no need to check ierr
         333             :
         334        1536 :    ierr = PIO_Inq_DimID(File, 'lev', lev_dimid)
         335        1536 :    ierr = PIO_Inq_dimlen(File, lev_dimid, flev)
         336        1536 :    if (plev /= flev) then
         337           0 :       write(iulog,*) routine//': ERROR: file lev does not match model. lev (file, model):',flev, plev
         338           0 :       call endrun(routine//': ERROR: file lev does not match model.')
         339             :    end if
         340             :
         341        1536 :    ierr = PIO_Inq_DimID(File, 'ilev', lev_dimid)
         342        1536 :    ierr = PIO_Inq_dimlen(File, lev_dimid, filev)
         343        1536 :    if (plevp /= filev) then
         344           0 :       write(iulog,*) routine//':ERROR: file ilev does not match model plevp (file, model):',filev, plevp
         345           0 :       call endrun(routine//':ERROR: file ilev does not match model.')
         346             :    end if
         347             :
         348        1536 :    ierr = pio_inq_varid(File, 'hyai', hyai_desc)
         349        1536 :    ierr = pio_inq_varid(File, 'hyam', hyam_desc)
         350        1536 :    ierr = pio_inq_varid(File, 'hybi', hybi_desc)
         351        1536 :    ierr = pio_inq_varid(File, 'hybm', hybm_desc)
         352             :
         353        1536 :    ierr = pio_get_var(File, hyai_desc, hyai)
         354        1536 :    ierr = pio_get_var(File, hybi_desc, hybi)
         355        1536 :    ierr = pio_get_var(File, hyam_desc, hyam)
         356        1536 :    ierr = pio_get_var(File, hybm_desc, hybm)
         357             :
         358        1536 :    if (masterproc) then
         359           2 :       write(iulog,*) routine//': read hyai, hybi, hyam, hybm'
         360             :    end if
         361             :
         362             :    ! Check whether file contains value for P0.  If it does then use it
         363             :
         364             :    ! Set PIO to return error codes.
         365        1536 :    call pio_seterrorhandling(file, PIO_BCAST_ERROR, pio_errtype)
         366             :
         367        1536 :    ierr = pio_inq_varid(file, 'P0', p0_desc)
         368        1536 :    if (ierr == PIO_NOERR) then
         369         768 :       ierr = pio_get_var(file, p0_desc, ps0)
         370         768 :       if (ierr /= PIO_NOERR) then
         371           0 :          call endrun(routine//': reading P0.')
         372             :       end if
         373         768 :       psr = ps0
         374             :
         375         768 :       if (masterproc) then
         376           1 :          write(iulog,*) routine//': read P0 value: ', ps0
         377             :       end if
         378             :
         379             :    end if
         380             :
         381             :    ! Put the error handling back the way it was
         382        1536 :    call pio_seterrorhandling(file, pio_errtype)
         383             :
         384             : #if ( defined OFFLINE_DYN )
         385             :    ! make sure top interface is non zero for fv dycore
         386             :    if (hyai(1) .eq. 0._r8) then
         387             :       if (hybm(1) .ne. 0.0_r8) then
         388             :          hyai(1) = hybm(1)*1.e-2_r8
         389             :       else if (hyam(1) .ne. 0.0_r8) then
         390             :          hyai(1) = hyam(1)*1.e-2_r8
         391             :       else
         392             :          call endrun('Not able to set hyai(1) to non-zero.')
         393             :       end if
         394             :    end if
         395             : #endif
         396             :
         397        1536 : end subroutine hycoef_read
         398             :
         399             : !=======================================================================
         400             :
         401             : end module hycoef

Generated by: LCOV version 1.14