LCOV - code coverage report
Current view: | top level - physics/cam - phys_grid_ctem.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 65 | 178 | 36.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 6 | 7 | 85.7 %

              Line data    Source code
           1             : !----------------------------------------------------------------------------------
           2             : ! circulation diagnostics -- terms of the Transformed Eulerian Mean (TEM) equation
           3             : !
           4             : !----------------------------------------------------------------------------------
           5             : module phys_grid_ctem
           6             :   use shr_kind_mod,  only: r8 => shr_kind_r8
           7             :   use ppgrid,        only: begchunk, endchunk, pcols, pver
           8             :   use physics_types, only: physics_state
           9             :   use cam_history,   only: addfld, outfld
          10             :   use zonal_mean_mod,only: ZonalAverage_t, ZonalMean_t
          11             :   use physconst,     only: pi
          12             :   use cam_logfile,   only: iulog
          13             :   use cam_abortutils,only: endrun, handle_allocate_error
          14             :   use namelist_utils,only: find_group_name
          15             :   use spmd_utils,    only: masterproc, mpi_integer, masterprocid, mpicom
          16             :   use time_manager,  only: get_step_size, get_nstep
          17             :
          18             :   use shr_const_mod, only: rgas => shr_const_rgas ! J/K/kmole
          19             :   use shr_const_mod, only: grav => shr_const_g ! m/s2
          20             :   use air_composition, only: mbarv ! g/mole
          21             :   use string_utils,  only: int2str
          22             :
          23             :   implicit none
          24             :
          25             :   private
          26             :   public :: phys_grid_ctem_readnl
          27             :   public :: phys_grid_ctem_reg
          28             :   public :: phys_grid_ctem_init
          29             :   public :: phys_grid_ctem_diags
          30             :   public :: phys_grid_ctem_final
          31             :
          32             :   type(ZonalMean_t) :: ZMobj
          33             :   type(ZonalAverage_t) :: ZAobj
          34             :
          35             :   integer :: nzalat = -huge(1)
          36             :   integer :: nzmbas = -huge(1)
          37             :
          38             :   integer :: ntimesteps = -huge(1) ! number of time steps bewteen TEM calculations
          39             :
          40             :   logical :: do_tem_diags = .false.
          41             :
          42             : contains
          43             :
          44             :   !------------------------------------------------------------------------------
          45             :   !------------------------------------------------------------------------------
          46        1536 :   subroutine phys_grid_ctem_readnl(nlfile)
          47             :     character(len=*), intent(in) :: nlfile
          48             :     integer :: ierr, unitn
          49             :
          50             :     character(len=*), parameter :: prefix = 'phys_grid_ctem_readnl: '
          51             :     real(r8) :: dtime
          52             :
          53             :     integer :: phys_grid_ctem_zm_nbas
          54             :     integer :: phys_grid_ctem_za_nlat
          55             :     integer :: phys_grid_ctem_nfreq
          56             :
          57             :     namelist /phys_grid_ctem_opts/ phys_grid_ctem_zm_nbas, phys_grid_ctem_za_nlat, phys_grid_ctem_nfreq
          58             :
          59        1536 :     phys_grid_ctem_zm_nbas = 0
          60        1536 :     phys_grid_ctem_za_nlat = 0
          61        1536 :     phys_grid_ctem_nfreq = 0
          62             :
          63             :     ! Read in namelist values
          64             :     !------------------------
          65        1536 :     if(masterproc) then
          66           2 :        open(newunit=unitn, file=trim(nlfile), status='old')
          67           2 :        call find_group_name(unitn, 'phys_grid_ctem_opts', status=ierr)
          68           2 :        if(ierr == 0) then
          69           0 :           read(unitn,phys_grid_ctem_opts,iostat=ierr)
          70           0 :           if(ierr /= 0) then
          71           0 :              call endrun(prefix//'ERROR reading namelist')
          72             :           end if
          73             :        end if
          74           2 :        close(unitn)
          75             :     end if
          76             :
          77        1536 :     call MPI_bcast(phys_grid_ctem_zm_nbas, 1, mpi_integer, masterprocid, mpicom, ierr)
          78        1536 :     if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_zm_nbas')
          79        1536 :     call MPI_bcast(phys_grid_ctem_za_nlat, 1, mpi_integer, masterprocid, mpicom, ierr)
          80        1536 :     if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_za_nlat')
          81        1536 :     call MPI_bcast(phys_grid_ctem_nfreq,   1, mpi_integer, masterprocid, mpicom, ierr)
          82        1536 :     if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_nfreq')
          83             :
          84        1536 :     do_tem_diags = .false.
          85        1536 :     if (phys_grid_ctem_nfreq/=0) then
          86           0 :        if (.not.(phys_grid_ctem_zm_nbas>0 .and. phys_grid_ctem_za_nlat>0)) then
          87             :           call endrun(prefix//'inconsistent phys_grid_ctem namelist settings -- phys_grid_ctem_zm_nbas=' &
          88           0 :                       //int2str(phys_grid_ctem_zm_nbas)//', phys_grid_ctem_za_nlat='//int2str(phys_grid_ctem_za_nlat))
          89             :        end if
          90           0 :        if (phys_grid_ctem_nfreq>0) then
          91           0 :           ntimesteps = phys_grid_ctem_nfreq
          92             :        else
          93           0 :           dtime = get_step_size()
          94           0 :           ntimesteps = nint( -phys_grid_ctem_nfreq*3600._r8/dtime )
          95             :        end if
          96           0 :        if (ntimesteps<1) then
          97             :           call endrun(prefix//'invalid ntimesteps -- phys_grid_ctem_nfreq needs to be a larger negative value ' &
          98           0 :                             //'or the model time step needs to be shorter')
          99             :        end if
         100           0 :        do_tem_diags = .true.
         101             :     end if
         102             :
         103        1536 :     if (masterproc) then
         104           2 :        if (do_tem_diags) then
         105           0 :           write(iulog,*) 'TEM diagnostics will be calculated every ',ntimesteps,' time steps'
         106           0 :           write(iulog,*) ' phys_grid_ctem_zm_nbas = ', phys_grid_ctem_zm_nbas
         107           0 :           write(iulog,*) ' phys_grid_ctem_za_nlat = ', phys_grid_ctem_za_nlat
         108           0 :           write(iulog,*) ' phys_grid_ctem_nfreq = ', phys_grid_ctem_nfreq
         109             :        else
         110           2 :           write(iulog,*) 'TEM diagnostics will not be performed'
         111             :        end if
         112             :     endif
         113             :
         114        1536 :     if (do_tem_diags) then
         115           0 :        nzalat = phys_grid_ctem_za_nlat
         116           0 :        nzmbas = phys_grid_ctem_zm_nbas
         117             :     end if
         118             :
         119        1536 :   end subroutine phys_grid_ctem_readnl
         120             :
         121             :   !-----------------------------------------------------------------------------
         122             :   !-----------------------------------------------------------------------------
         123        1536 :   subroutine phys_grid_ctem_reg
         124             :
         125             :     use cam_grid_support, only: horiz_coord_t, horiz_coord_create, iMap, cam_grid_register
         126             :
         127             :     type(horiz_coord_t), pointer :: zalon_coord
         128             :     type(horiz_coord_t), pointer :: zalat_coord
         129        1536 :     integer(iMap),       pointer :: grid_map(:,:)
         130             :
         131        3072 :     real(r8) :: zalats(nzalat)
         132        3072 :     real(r8) :: area(nzalat)
         133             :     real(r8) :: zalons(1)
         134             :     real(r8) :: dlatrad, dlatdeg, lat1, lat2
         135             :     real(r8) :: total_area
         136             :     real(r8) :: total_wght
         137             :     integer :: j, astat
         138             :
         139             :     real(r8), parameter :: latdeg0 = -90._r8
         140             :     real(r8), parameter :: latrad0 = -pi*0.5_r8
         141             :     real(r8), parameter :: fourpi = pi*4._r8
         142             :
         143             :     integer, parameter :: ctem_zavg_phys_decomp = 333 ! Must be unique within CAM
         144             :
         145        1536 :     if (.not.do_tem_diags) return
         146             :
         147           0 :     nullify(zalat_coord)
         148           0 :     nullify(zalon_coord)
         149           0 :     nullify(grid_map)
         150             :
         151           0 :     zalons(1) = 0._r8
         152             :
         153           0 :     dlatrad = pi/real(nzalat,kind=r8)
         154           0 :     dlatdeg = 180._r8/real(nzalat,kind=r8)
         155           0 :     total_area = 0._r8
         156           0 :     total_wght = 0._r8
         157             :
         158             :     ! calculate latitudes and areas of zonal average grid boxes
         159           0 :     do j = 1,nzalat
         160           0 :        zalats(j) = latdeg0 + (real(j,kind=r8)-0.5_r8)*dlatdeg
         161           0 :        lat1 = latrad0 + real(j-1,kind=r8)*dlatrad
         162           0 :        lat2 = latrad0 + real(j  ,kind=r8)*dlatrad
         163           0 :        area(j) = 2._r8*pi*(sin(lat2)-sin(lat1))
         164           0 :        total_area = total_area + area(j)
         165           0 :        total_wght = total_wght + 0.5_r8*(sin(lat2)-sin(lat1))
         166             :     end do
         167             :
         168             :     ! sanity check
         169           0 :     if ( abs(1._r8-total_wght)>1.e-12_r8 .or. abs(fourpi-total_area)>1.e-12_r8 ) then
         170           0 :        call endrun('phys_grid_ctem_reg: problem with area/wght calc')
         171             :     end if
         172             :
         173             :     ! initialize zonal-average and zonal-mean utility objects
         174           0 :     call ZAobj%init(zalats,area,nzalat,GEN_GAUSSLATS=.false.)
         175           0 :     call ZMobj%init(nzmbas)
         176             :
         177             :     ! Zonal average grid for history fields
         178             :
         179           0 :     zalat_coord => horiz_coord_create('zalat', '', nzalat, 'latitude', 'degrees_north', 1, nzalat, zalats)
         180           0 :     zalon_coord => horiz_coord_create('zalon', '', 1, 'longitude', 'degrees_east', 1, 1, zalons)
         181             :
         182             :     ! grid decomposition map
         183           0 :     allocate(grid_map(4,nzalat), stat=astat)
         184           0 :     call handle_allocate_error(astat, 'phys_grid_ctem_reg', 'grid_map')
         185             :
         186           0 :     do j = 1,nzalat
         187           0 :        grid_map(1,j) = 1
         188           0 :        grid_map(2,j) = j
         189           0 :        if (masterproc) then
         190           0 :           grid_map(3,j) = 1
         191           0 :           grid_map(4,j) = j
         192             :        else
         193           0 :           grid_map(3,j) = 0
         194           0 :           grid_map(4,j) = 0
         195             :        end if
         196             :     end do
         197             :
         198             :     ! register the zonal average grid
         199             :     call cam_grid_register('ctem_zavg_phys', ctem_zavg_phys_decomp, zalat_coord, zalon_coord, grid_map, &
         200           0 :                            unstruct=.false., zonal_grid=.true.)
         201             :
         202        1536 :   end subroutine phys_grid_ctem_reg
         203             :
         204             :   !-----------------------------------------------------------------------------
         205             :   !-----------------------------------------------------------------------------
         206        1536 :   subroutine phys_grid_ctem_init
         207             :
         208        1536 :     if (.not.do_tem_diags) return
         209             :
         210           0 :     call addfld ('Uzm',  (/'lev'/), 'A','m s-1',  'Zonal-Mean zonal wind', gridname='ctem_zavg_phys' )
         211           0 :     call addfld ('Vzm',  (/'lev'/), 'A','m s-1',  'Zonal-Mean meridional wind', gridname='ctem_zavg_phys' )
         212           0 :     call addfld ('Wzm',  (/'lev'/), 'A','m s-1',  'Zonal-Mean vertical wind', gridname='ctem_zavg_phys' )
         213           0 :     call addfld ('THzm', (/'lev'/), 'A','K',      'Zonal-Mean potential temp', gridname='ctem_zavg_phys' )
         214           0 :     call addfld ('VTHzm',(/'lev'/), 'A','K m s-1','Meridional Heat Flux:', gridname='ctem_zavg_phys')
         215           0 :     call addfld ('WTHzm',(/'lev'/), 'A','K m s-1','Vertical Heat Flux:', gridname='ctem_zavg_phys')
         216           0 :     call addfld ('UVzm', (/'lev'/), 'A','m2 s-2', 'Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys')
         217           0 :     call addfld ('UWzm', (/'lev'/), 'A','m2 s-2', 'Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys')
         218           0 :     call addfld ('THphys',(/'lev'/), 'A', 'K',    'Potential temp', gridname='physgrid' )
         219             :
         220        1536 :   end subroutine phys_grid_ctem_init
         221             :
         222             :   !-----------------------------------------------------------------------------
         223             :   !-----------------------------------------------------------------------------
         224      370944 :   subroutine phys_grid_ctem_diags(phys_state)
         225             :     type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
         226             :
         227             :     character(len=*), parameter :: prefix = 'phys_grid_ctem_diags: '
         228             :
         229      741888 :     real(r8) :: u(pcols,pver,begchunk:endchunk)
         230      741888 :     real(r8) :: v(pcols,pver,begchunk:endchunk)
         231      741888 :     real(r8) :: w(pcols,pver,begchunk:endchunk)
         232             :
         233      741888 :     real(r8) :: uzm(pcols,pver,begchunk:endchunk)
         234      741888 :     real(r8) :: vzm(pcols,pver,begchunk:endchunk)
         235      741888 :     real(r8) :: wzm(pcols,pver,begchunk:endchunk)
         236             :
         237      741888 :     real(r8) :: ud(pcols,pver,begchunk:endchunk)
         238      741888 :     real(r8) :: vd(pcols,pver,begchunk:endchunk)
         239      741888 :     real(r8) :: wd(pcols,pver,begchunk:endchunk)
         240      741888 :     real(r8) :: thd(pcols,pver,begchunk:endchunk)
         241             :
         242      741888 :     real(r8) :: uvp(pcols,pver,begchunk:endchunk)
         243      741888 :     real(r8) :: uwp(pcols,pver,begchunk:endchunk)
         244      741888 :     real(r8) :: vthp(pcols,pver,begchunk:endchunk)
         245      741888 :     real(r8) :: wthp(pcols,pver,begchunk:endchunk)
         246             :
         247             :     integer  :: lchnk, ncol, j, k
         248             :
         249             :     ! potential temperature
         250      741888 :     real(r8) :: theta(pcols,pver,begchunk:endchunk)
         251      741888 :     real(r8) :: thzm(pcols,pver,begchunk:endchunk)
         252             :
         253      741888 :     real(r8) :: uvza(nzalat,pver)
         254      741888 :     real(r8) :: uwza(nzalat,pver)
         255      741888 :     real(r8) :: vthza(nzalat,pver)
         256      741888 :     real(r8) :: wthza(nzalat,pver)
         257             :
         258      741888 :     real(r8) :: uza(nzalat,pver)
         259      741888 :     real(r8) :: vza(nzalat,pver)
         260      741888 :     real(r8) :: wza(nzalat,pver)
         261      370944 :     real(r8) :: thza(nzalat,pver)
         262             :
         263             :     real(r8) :: sheight(pcols,pver) ! pressure scale height (m)
         264             :
         265      370944 :     if (.not.do_calc()) return
         266             :
         267           0 :     do lchnk = begchunk,endchunk
         268             :
         269           0 :        ncol = phys_state(lchnk)%ncol
         270             :
         271             :        ! scale height
         272           0 :        sheight(:ncol,:) = phys_state(lchnk)%t(:ncol,:) * rgas / ( mbarv(:ncol,:,lchnk) * grav ) ! meters
         273             :
         274             :        ! potential temperature
         275           0 :        theta(:ncol,:,lchnk) = phys_state(lchnk)%t(:ncol,:) * phys_state(lchnk)%exner(:ncol,:)
         276             :
         277             :        ! vertical velocity
         278           0 :        w(:ncol,:,lchnk) = -sheight(:ncol,:) *  phys_state(lchnk)%omega(:ncol,:) / phys_state(lchnk)%pmid(:ncol,:)
         279             :
         280           0 :        u(:ncol,:,lchnk) =  phys_state(lchnk)%u(:ncol,:)
         281           0 :        v(:ncol,:,lchnk) =  phys_state(lchnk)%v(:ncol,:)
         282             :
         283             :     end do
         284             :
         285             :     ! zonal means evaluated on the physics grid (3D) to be used in the deviations calculation below
         286           0 :     uzm(:,:,:) = zmean_fld(u(:,:,:))
         287           0 :     vzm(:,:,:) = zmean_fld(v(:,:,:))
         288           0 :     wzm(:,:,:) = zmean_fld(w(:,:,:))
         289           0 :     thzm(:,:,:) = zmean_fld(theta(:,:,:))
         290             :
         291             :     ! diagnostic output
         292           0 :     do lchnk = begchunk, endchunk
         293           0 :        call outfld( 'THphys', theta(:,:,lchnk), pcols, lchnk)
         294             :     end do
         295             :
         296           0 :     do lchnk = begchunk,endchunk
         297           0 :        ncol = phys_state(lchnk)%ncol
         298           0 :        do k = 1,pver
         299             :           ! zonal deviations
         300           0 :           thd(:ncol,k,lchnk) = theta(:ncol,k,lchnk) - thzm(:ncol,k,lchnk)
         301           0 :           ud(:ncol,k,lchnk) = u(:ncol,k,lchnk) - uzm(:ncol,k,lchnk)
         302           0 :           vd(:ncol,k,lchnk) = v(:ncol,k,lchnk) - vzm(:ncol,k,lchnk)
         303           0 :           wd(:ncol,k,lchnk) = w(:ncol,k,lchnk) - wzm(:ncol,k,lchnk)
         304             :           ! fluxes
         305           0 :           uvp(:ncol,k,lchnk) = ud(:ncol,k,lchnk) * vd(:ncol,k,lchnk)
         306           0 :           uwp(:ncol,k,lchnk) = ud(:ncol,k,lchnk) * wd(:ncol,k,lchnk)
         307           0 :           vthp(:ncol,k,lchnk) = vd(:ncol,k,lchnk) * thd(:ncol,k,lchnk)
         308           0 :           wthp(:ncol,k,lchnk) = wd(:ncol,k,lchnk) * thd(:ncol,k,lchnk)
         309             :        end do
         310             :     end do
         311             :
         312             :     ! evaluate and output fluxes on the zonal-average grid
         313           0 :     call ZAobj%binAvg(uvp, uvza)
         314           0 :     call ZAobj%binAvg(uwp, uwza)
         315           0 :     call ZAobj%binAvg(vthp, vthza)
         316           0 :     call ZAobj%binAvg(wthp, wthza)
         317             :
         318           0 :     if (any(abs(uvza)>1.e20_r8)) call endrun(prefix//'bad values in uvza')
         319           0 :     if (any(abs(uwza)>1.e20_r8)) call endrun(prefix//'bad values in uwza')
         320           0 :     if (any(abs(vthza)>1.e20_r8)) call endrun(prefix//'bad values in vthza')
         321           0 :     if (any(abs(wthza)>1.e20_r8)) call endrun(prefix//'bad values in wthza')
         322             :
         323           0 :     call ZAobj%binAvg(uzm, uza)
         324           0 :     call ZAobj%binAvg(vzm, vza)
         325           0 :     call ZAobj%binAvg(wzm, wza)
         326           0 :     call ZAobj%binAvg(thzm, thza)
         327             :
         328           0 :     if (any(abs(uza)>1.e20_r8)) call endrun(prefix//'bad values in uza')
         329           0 :     if (any(abs(vza)>1.e20_r8)) call endrun(prefix//'bad values in vza')
         330           0 :     if (any(abs(wza)>1.e20_r8)) call endrun(prefix//'bad values in wza')
         331           0 :     if (any(abs(thza)>1.e20_r8)) call endrun(prefix//'bad values in thza')
         332             :
         333             :     ! diagnostic output
         334           0 :     do j = 1,nzalat
         335           0 :        call outfld('Uzm',uza(j,:),1,j)
         336           0 :        call outfld('Vzm',vza(j,:),1,j)
         337           0 :        call outfld('Wzm',wza(j,:),1,j)
         338           0 :        call outfld('THzm',thza(j,:),1,j)
         339           0 :        call outfld('UVzm',uvza(j,:),1,j)
         340           0 :        call outfld('UWzm',uwza(j,:),1,j)
         341           0 :        call outfld('VTHzm',vthza(j,:),1,j)
         342           0 :        call outfld('WTHzm',wthza(j,:),1,j)
         343             :     end do
         344             :
         345             :   contains
         346             :
         347             :     !------------------------------------------------------------------------------
         348             :     ! utility function for evaluating 3D zonal mean fields
         349             :     !------------------------------------------------------------------------------
         350           0 :     function zmean_fld( fld ) result(fldzm)
         351             :
         352             :       real(r8), intent(in) :: fld(pcols,pver,begchunk:endchunk)
         353             :
         354             :       real(r8) :: fldzm(pcols,pver,begchunk:endchunk)
         355             :
         356           0 :       real(r8) :: Zonal_Bamp3d(nzmbas,pver)
         357             :
         358           0 :       call ZMobj%calc_amps(fld,Zonal_Bamp3d)
         359           0 :       call ZMobj%eval_grid(Zonal_Bamp3d,fldzm)
         360             :
         361           0 :     end function zmean_fld
         362             :
         363             :     !------------------------------------------------------------------------------
         364             :     ! utility function returns TRUE when time to update TEM diags
         365             :     !------------------------------------------------------------------------------
         366      370944 :     logical function do_calc()
         367             :
         368             :       integer :: nstep
         369      370944 :       nstep = get_nstep()
         370      370944 :       do_calc = do_tem_diags .and. mod(nstep,ntimesteps) == 0
         371             :
         372      370944 :     end function do_calc
         373             :
         374             :   end subroutine phys_grid_ctem_diags
         375             :
         376             :   !-----------------------------------------------------------------------------
         377             :   !-----------------------------------------------------------------------------
         378        1536 :   subroutine phys_grid_ctem_final
         379        1536 :     call ZAobj%final()
         380        1536 :     call ZMobj%final()
         381        1536 :   end subroutine phys_grid_ctem_final
         382             :
         383             : end module phys_grid_ctem

Generated by: LCOV version 1.14