module datm_shr_mod ! !USES: use shr_kind_mod , only : IN=>SHR_KIND_IN, R8=>SHR_KIND_R8 use shr_kind_mod , only : CS=>SHR_KIND_CS, CL=>SHR_KIND_CL use shr_const_mod , only : SHR_CONST_CDAY,SHR_CONST_TKFRZ,SHR_CONST_SPVAL use shr_file_mod , only : shr_file_getlogunit, shr_file_getunit, shr_file_freeunit use shr_sys_mod , only : shr_sys_flush, shr_sys_abort use shr_strdata_mod, only : shr_strdata_readnml use shr_dmodel_mod , only : shr_dmodel_mapset use seq_flds_mod , only : seq_flds_dom_coord, seq_flds_dom_other use shr_cal_mod , only : shr_cal_date2julian use shr_ncread_mod , only : shr_ncread_varExists, shr_ncread_varDimSizes, shr_ncread_field4dG use shr_strdata_mod, only : shr_strdata_type use mct_mod ! !PUBLIC TYPES: implicit none private ! except !-------------------------------------------------------------------------- ! Public interfaces !-------------------------------------------------------------------------- public :: datm_shr_getNextRadCDay public :: datm_shr_CORE2getFactors public :: datm_shr_TN460getFactors public :: datm_shr_eSat public :: datm_shr_read_namelists !-------------------------------------------------------------------------- ! Public data !-------------------------------------------------------------------------- ! input namelist variables character(CL) , public :: decomp ! decomp strategy character(CL) , public :: restfilm ! model restart file namelist character(CL) , public :: restfils ! stream restart file namelist character(CL) , public :: bias_correct ! true => send bias correction fields to coupler character(CL) , public :: anomaly_forcing(8) ! true => send anomaly forcing fields to coupler logical , public :: force_prognostic_true ! if true set prognostic true logical , public :: wiso_datm = .false. ! expect isotopic forcing from file? integer(IN) , public :: iradsw ! radiation interval character(CL) , public :: factorFn ! file containing correction factors logical , public :: presaero ! true => send valid prescribe aero fields to coupler ! variables obtained from namelist read character(CL) , public :: rest_file ! restart filename character(CL) , public :: rest_file_strm ! restart filename for streams character(CL) , public :: datamode ! mode character(len=*), public, parameter :: nullstr = 'undefined' !-------------------------------------------------------------------------- ! Private data !-------------------------------------------------------------------------- save !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CONTAINS !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine datm_shr_read_namelists(mpicom, my_task, master_task, & inst_index, inst_suffix, inst_name, & logunit, shrlogunit, SDATM, atm_present, atm_prognostic) ! !INPUT/OUTPUT PARAMETERS: integer(IN) , intent(in) :: mpicom ! mpi communicator integer(IN) , intent(in) :: my_task ! my task in mpi communicator mpicom integer(IN) , intent(in) :: master_task ! task number of master task integer , intent(in) :: inst_index ! number of current instance (ie. 1) character(len=16) , intent(in) :: inst_suffix ! char string associated with instance character(len=16) , intent(in) :: inst_name ! fullname of current instance (ie. "lnd_0001") integer(IN) , intent(in) :: logunit ! logging unit number integer(IN) , intent(in) :: shrlogunit ! original log unit and level type(shr_strdata_type) , intent(inout) :: SDATM logical , intent(out) :: atm_present ! flag logical , intent(out) :: atm_prognostic ! flag !--- local variables --- character(CL) :: fileName ! generic file name integer(IN) :: nunit ! unit number integer(IN) :: ierr ! error code !--- formats --- character(*), parameter :: F00 = "('(datm_comp_init) ',8a)" character(*), parameter :: F0L = "('(datm_comp_init) ',a, l2)" character(*), parameter :: F01 = "('(datm_comp_init) ',a,5i8)" character(*), parameter :: subName = "(shr_datm_read_namelists) " !------------------------------------------------------------------------------- !----- define namelist ----- namelist / datm_nml / & decomp, iradsw, factorFn, restfilm, restfils, presaero, bias_correct, & anomaly_forcing, force_prognostic_true, wiso_datm !---------------------------------------------------------------------------- ! Determine input filenamname !---------------------------------------------------------------------------- filename = "datm_in"//trim(inst_suffix) !---------------------------------------------------------------------------- ! Read datm_in !---------------------------------------------------------------------------- decomp = "1d" iradsw = 0 factorFn = 'null' restfilm = trim(nullstr) restfils = trim(nullstr) presaero = .false. force_prognostic_true = .false. if (my_task == master_task) then nunit = shr_file_getUnit() ! get unused unit number open (nunit,file=trim(filename),status="old",action="read") read (nunit,nml=datm_nml,iostat=ierr) close(nunit) call shr_file_freeUnit(nunit) if (ierr > 0) then write(logunit,F01) 'ERROR: reading input namelist, '//trim(filename)//' iostat=',ierr call shr_sys_abort(subName//': namelist read error '//trim(filename)) end if write(logunit,F00)' decomp = ',trim(decomp) write(logunit,F01)' iradsw = ',iradsw write(logunit,F00)' factorFn = ',trim(factorFn) write(logunit,F00)' restfilm = ',trim(restfilm) write(logunit,F00)' restfils = ',trim(restfils) write(logunit,F0L)' presaero = ',presaero write(logunit,F0L)' force_prognostic_true = ',force_prognostic_true write(logunit,F0L)' wiso_datm = ', wiso_datm write(logunit,F01) 'inst_index = ',inst_index write(logunit,F00) 'inst_name = ',trim(inst_name) write(logunit,F00) 'inst_suffix = ',trim(inst_suffix) call shr_sys_flush(logunit) endif call shr_mpi_bcast(decomp,mpicom,'decomp') call shr_mpi_bcast(iradsw,mpicom,'iradsw') call shr_mpi_bcast(factorFn,mpicom,'factorFn') call shr_mpi_bcast(restfilm,mpicom,'restfilm') call shr_mpi_bcast(restfils,mpicom,'restfils') call shr_mpi_bcast(presaero,mpicom,'presaero') call shr_mpi_bcast(force_prognostic_true,mpicom,'force_prognostic_true') call shr_mpi_bcast(wiso_datm, mpicom, 'wiso_datm') rest_file = trim(restfilm) rest_file_strm = trim(restfils) !---------------------------------------------------------------------------- ! Read dshr namelist !---------------------------------------------------------------------------- call shr_strdata_readnml(SDATM, trim(filename), mpicom=mpicom) call shr_sys_flush(shrlogunit) ! Validate mode datamode = trim(SDATM%dataMode) if (trim(datamode) == 'NULL' .or. & trim(datamode) == 'CORE2_NYF' .or. & trim(datamode) == 'CORE2_IAF' .or. & trim(datamode) == 'CORE_IAF_JRA' .or. & trim(datamode) == 'CLMNCEP' .or. & trim(datamode) == 'COPYALL' ) then if (my_task == master_task) then write(logunit,F00) ' datm datamode = ',trim(datamode) call shr_sys_flush(logunit) end if else write(logunit,F00) ' ERROR illegal datm datamode = ',trim(datamode) call shr_sys_abort() endif !---------------------------------------------------------------------------- ! Determine present and prognostic flag !---------------------------------------------------------------------------- atm_present = .false. atm_prognostic = .false. if (force_prognostic_true) then atm_present = .true. atm_prognostic = .true. endif if (trim(datamode) /= 'NULL') then atm_present = .true. end if end subroutine datm_shr_read_namelists !=============================================================================== real(R8) function datm_shr_getNextRadCDay( ymd, tod, stepno, dtime, iradsw, calendar ) ! !DESCRIPTION: ! Return the calendar day of the next radiation time-step. ! General Usage: nextswday = datm_shr_getNextRadCDay(curr_date) implicit none ! !INPUT/OUTPUT PARAMETERS: integer(IN), intent(IN) :: ymd integer(IN), intent(IN) :: tod integer(IN), intent(IN) :: stepno integer(IN), intent(IN) :: dtime integer(IN), intent(IN) :: iradsw character(*),intent(in) :: calendar !----- local ----- real(R8) :: nextsw_cday real(R8) :: julday integer :: liradsw character(*),parameter :: subName = '(datm_shr_getNextRadCDay) ' !------------------------------------------------------------------------------- liradsw = iradsw if (liradsw < 0) liradsw = nint((-liradsw *3600._r8)/dtime) call shr_cal_date2julian(ymd,tod,julday,calendar) if (liradsw > 1) then if (mod(stepno+1,liradsw) == 0 .and. stepno > 0) then nextsw_cday = julday + 2*dtime/SHR_CONST_CDAY else nextsw_cday = -1._r8 end if else nextsw_cday = julday + dtime/SHR_CONST_CDAY end if datm_shr_getNextRadCDay = nextsw_cday end function datm_shr_getNextRadCDay !=============================================================================== subroutine datm_shr_CORE2getFactors(fileName,windF,winddF,qsatF,mpicom,compid, & gsmap,ggrid,nxg,nyg) implicit none !--- arguments --- character(*) ,intent(in) :: fileName ! file name string real(R8) ,intent(inout) :: windF(:) ! wind adjustment factor real(R8) ,intent(inout) :: winddF(:) ! wind adjustment factor real(R8) ,intent(inout) :: qsatF(:) ! rel humidty adjustment factors integer(IN) ,intent(in) :: mpicom ! mpi comm integer(IN) ,intent(in) :: compid ! mct compid type(mct_gsmap) ,intent(in) :: gsmap ! decomp of wind,windd,qsat type(mct_ggrid) ,intent(in) :: ggrid ! ggrid of grid info integer(IN) ,intent(in) :: nxg ! size of input grid integer(IN) ,intent(in) :: nyg ! size of input grid !--- local --- integer(IN) :: my_task,logunit,ier !--- formats --- character(*),parameter :: subName = '(datm_shr_CORE2getFactors) ' character(*),parameter :: F00 = "('(datm_shr_CORE2getFactors) ',4a) " !------------------------------------------------------------------------------- call MPI_COMM_RANK(mpicom,my_task,ier) call shr_file_getLogUnit(logunit) if (my_task == 0) then !--- verify necessary data is in input file --- if ( .not. shr_ncread_varExists(fileName,'lat' ) & .or. .not. shr_ncread_varExists(fileName,'lon' ) & .or. .not. shr_ncread_varExists(fileName,'mask' ) & .or. .not. shr_ncread_varExists(fileName,'windFactor') & .or. .not. shr_ncread_varExists(fileName,'winddFactor') & .or. .not. shr_ncread_varExists(fileName,'qsatFactor') ) then write(logunit,F00) "ERROR: invalid correction factor data file" call shr_sys_abort(subName//"invalid correction factor data file") end if endif call datm_shr_getFactors(fileName,windF,winddF,qsatF,mpicom,compid, & gsmap,ggrid,nxg,nyg) end subroutine datm_shr_CORE2getFactors !=============================================================================== subroutine datm_shr_TN460getFactors(fileName,windF,qsatF,mpicom,compid, & gsmap,ggrid,nxg,nyg) implicit none !--- arguments --- character(*) ,intent(in) :: fileName ! file name string real(R8) ,intent(inout) :: windF(:) ! wind adjustment factor real(R8) ,intent(inout) :: qsatF(:) ! rel humidty adjustment factors integer(IN) ,intent(in) :: mpicom ! mpi comm integer(IN) ,intent(in) :: compid ! mct compid type(mct_gsmap) ,intent(in) :: gsmap ! decomp of wind,windd,qsat type(mct_ggrid) ,intent(in) :: ggrid ! ggrid of grid info integer(IN) ,intent(in) :: nxg ! size of input grid integer(IN) ,intent(in) :: nyg ! size of input grid !--- local --- integer(IN) :: my_task,logunit,ier real(R8),pointer :: winddF(:) ! wind adjustment factor !--- formats --- character(*),parameter :: subName = '(datm_shr_TN460getFactors) ' character(*),parameter :: F00 = "('(datm_shr_TN460getFactors) ',4a) " !------------------------------------------------------------------------------- call MPI_COMM_RANK(mpicom,my_task,ier) call shr_file_getLogUnit(logunit) if (my_task == 0) then !--- verify necessary data is in input file --- if ( .not. shr_ncread_varExists(fileName,'lat' ) & .or. .not. shr_ncread_varExists(fileName,'lon' ) & .or. .not. shr_ncread_varExists(fileName,'mask' ) & .or. .not. shr_ncread_varExists(fileName,'windFactor') & .or. .not. shr_ncread_varExists(fileName,'qsatFactor') ) then write(logunit,F00) "ERROR: invalid correction factor data file" call shr_sys_abort(subName//"invalid correction factor data file") end if endif call datm_shr_getFactors(fileName,windF,winddF,qsatF,mpicom,compid, & gsmap,ggrid,nxg,nyg) end subroutine datm_shr_TN460getFactors !=============================================================================== subroutine datm_shr_getFactors(fileName,windF,winddF,qsatF,mpicom,compid, & gsmapo,ggrido,nxgo,nygo) use shr_map_mod implicit none !--- arguments --- character(*) ,intent(in) :: fileName ! file name string real(R8) ,intent(inout) :: windF(:) ! wind adjustment factor real(R8) ,intent(inout) :: winddF(:) ! wind adjustment factor real(R8) ,intent(inout) :: qsatF(:) ! rel humidty adjustment factors integer(IN) ,intent(in) :: mpicom ! mpi comm integer(IN) ,intent(in) :: compid ! mct compid type(mct_gsmap) ,intent(in) :: gsmapo ! decomp of wind,windd,qsat type(mct_ggrid) ,intent(in) :: ggrido ! ggrid of grid info integer(IN) ,intent(in) :: nxgo ! size of input grid integer(IN) ,intent(in) :: nygo ! size of input grid !--- data that describes the local model domain --- integer(IN) :: ni0,nj0 ! dimensions of global bundle0 integer(IN) :: i,j,n ! generic indicies integer(IN) :: my_task ! local pe number integer(IN) :: ier ! error code integer(IN) :: logunit ! logging unit type(mct_ggrid) :: ggridi ! input file grid type(mct_ggrid) :: ggridoG ! output grid gathered type(mct_gsmap) :: gsmapi ! input file gsmap type(mct_sMatp) :: smatp ! sparse matrix weights type(mct_avect) :: avi ! input attr vect type(mct_avect) :: avo ! output attr vect integer(IN) :: lsizei ! local size of input integer(IN) :: lsizeo ! local size of output integer(IN),pointer :: start(:) ! start list integer(IN),pointer :: length(:) ! length list integer(IN) :: gsizei ! input global size integer(IN) :: numel ! number of elements in start list real(R8) :: dadd ! lon correction logical :: domap ! map or not integer(IN) :: klon,klat ! lon lat fld index !--- temp arrays for data input --- real(R8) ,allocatable :: tempR4D(:,:,:,:) ! 4D data array real(R8) ,pointer :: tempR1D(:) ! 1D data array integer(IN),allocatable :: tempI4D(:,:,:,:) ! 4D data array !--- formats --- character(*),parameter :: subName = '(datm_shr_getFactors) ' character(*),parameter :: F00 = "('(datm_shr_getFactors) ',4a) " character(*),parameter :: F01 = "('(datm_shr_getFactors) ',a,2i5)" character(*),parameter :: F02 = "('(datm_shr_getFactors) ',a,6e12.3)" !------------------------------------------------------------------------------- ! Note: gsmapi is all gridcells on root pe !------------------------------------------------------------------------------- call MPI_COMM_RANK(mpicom,my_task,ier) call shr_file_getLogUnit(logunit) ni0 = 0 nj0 = 0 allocate(start(1),length(1)) start = 0 length = 0 numel = 0 !---------------------------------------------------------------------------- ! read in and map global correction factors !---------------------------------------------------------------------------- if (my_task == 0) then !--- verify necessary data is in input file --- if ( .not. shr_ncread_varExists(fileName,'lat' ) & .or. .not. shr_ncread_varExists(fileName,'lon' ) & .or. .not. shr_ncread_varExists(fileName,'mask' ) & .or. .not. shr_ncread_varExists(fileName,'windFactor') ) then write(logunit,F00) "ERROR: invalid correction factor data file" call shr_sys_abort(subName//"invalid correction factor data file") end if call shr_ncread_varDimSizes(fileName,"windFactor",ni0,nj0) start = 1 length = ni0*nj0 numel = 1 endif call shr_mpi_bcast(ni0,mpicom,subname//' ni0') call shr_mpi_bcast(nj0,mpicom,subname//' nj0') gsizei = ni0*nj0 !--- allocate datatypes for input data --- call mct_gsmap_init(gsmapi,start,length,0,mpicom,compid,gsize=gsizei,numel=numel) deallocate(start,length) lsizei = mct_gsmap_lsize(gsmapi,mpicom) lsizeo = mct_gsmap_lsize(gsmapo,mpicom) call mct_gGrid_init(GGrid=gGridi, CoordChars=trim(seq_flds_dom_coord), & OtherChars=trim(seq_flds_dom_other), lsize=lsizei ) call mct_aVect_init(avi,rList="wind:windd:qsat",lsize=lsizei) avi%rAttr = SHR_CONST_SPVAL !--- gather output grid for map logic --- call mct_ggrid_gather(ggrido, ggridoG, gsmapo, 0, mpicom) if (my_task == 0) then allocate(tempR1D(ni0*nj0)) !--- read domain data: lon --- allocate(tempR4D(ni0,1,1,1)) call shr_ncread_field4dG(fileName,'lon' ,rfld=tempR4D) !--- needs to be monotonically increasing, add 360 at wraparound+ --- dadd = 0.0_R8 do i = 2,ni0 if (tempR4D(i-1,1,1,1) > tempR4D(i,1,1,1)) dadd = 360.0_R8 tempR4D(i,1,1,1) = tempR4D(i,1,1,1) + dadd enddo n = 0 do j=1,nj0 do i=1,ni0 n = n + 1 tempR1D(n) = tempR4D(i,1,1,1) end do end do deallocate(tempR4D) call mct_gGrid_importRattr(gGridi,'lon',tempR1D,lsizei) !--- read domain data: lat --- allocate(tempR4D(nj0,1,1,1)) call shr_ncread_field4dG(fileName,'lat' ,rfld=tempR4D) n = 0 do j=1,nj0 do i=1,ni0 n = n + 1 tempR1D(n) = tempR4D(j,1,1,1) end do end do deallocate(tempR4D) call mct_gGrid_importRattr(gGridi,'lat',tempR1D,lsizei) !--- read domain mask--- allocate(tempI4D(ni0,nj0,1,1)) call shr_ncread_field4dG(fileName,'mask',ifld=tempI4D) n = 0 do j=1,nj0 do i=1,ni0 n = n + 1 tempR1D(n) = real(tempI4D(i,j,1,1),R8) end do end do deallocate(tempI4D) call mct_gGrid_importRattr(gGridi,'mask',tempR1D,lsizei) !--- read bundle data: wind factor --- allocate(tempR4D(ni0,nj0,1,1)) if (shr_ncread_varExists(fileName,'windFactor') ) then call shr_ncread_field4dG(fileName,'windFactor',rfld=tempR4D) n = 0 do j=1,nj0 do i=1,ni0 n = n + 1 tempR1D(n) = tempR4D(i,j,1,1) end do end do call mct_aVect_importRattr(avi,'wind',tempR1D,lsizei) endif !--- read bundle data: windd factor --- if (shr_ncread_varExists(fileName,'winddFactor') ) then call shr_ncread_field4dG(fileName,'winddFactor',rfld=tempR4D) n = 0 do j=1,nj0 do i=1,ni0 n = n + 1 tempR1D(n) = tempR4D(i,j,1,1) end do end do call mct_aVect_importRattr(avi,'windd',tempR1D,lsizei) endif !--- read bundle data: qsat factor --- if (shr_ncread_varExists(fileName,'qsatFactor') ) then call shr_ncread_field4dG(fileName,'qsatFactor',rfld=tempR4D) n = 0 do j=1,nj0 do i=1,ni0 n = n + 1 tempR1D(n) = tempR4D(i,j,1,1) end do end do call mct_aVect_importRattr(avi,'qsat',tempR1D,lsizei) endif deallocate(tempR4D) deallocate(tempR1D) domap = .false. if (ni0 /= nxgo .or. nj0 /= nygo) then domap = .true. else klon = mct_aVect_indexRA(ggridi%data,'lon') klat = mct_aVect_indexRA(ggrido%data,'lat') do n = 1,lsizei if (abs(ggridi%data%rAttr(klon,n)-ggridoG%data%rAttr(klon,n)) > 0.01_R8) domap=.true. if (abs(ggridi%data%rAttr(klat,n)-ggridoG%data%rAttr(klat,n)) > 0.01_R8) domap=.true. enddo endif call mct_gGrid_clean(ggridoG) endif call shr_mpi_bcast(domap,mpicom,subname//' domap') if (domap) then call shr_dmodel_mapSet(smatp,ggridi,gsmapi,ni0 ,nj0 , & ggrido,gsmapo,nxgo,nygo, & 'datmfactor',shr_map_fs_remap,shr_map_fs_bilinear,shr_map_fs_srcmask, & shr_map_fs_scalar,compid,mpicom,'Xonly') call mct_aVect_init(avo,avi,lsizeo) call mct_sMat_avMult(avi,smatp,avo) call mct_sMatP_clean(smatp) else call mct_aVect_scatter(avi,avo,gsmapo,0,mpicom) endif !--- fill the interface arrays, only if they are the right size --- allocate(tempR1D(lsizeo)) if (size(windF ) >= lsizeo) then call mct_aVect_exportRattr(avo,'wind' ,tempR1D,lsizeo) windF = tempR1D endif if (size(winddF) >= lsizeo) then call mct_aVect_exportRattr(avo,'windd',tempR1D,lsizeo) winddF = tempR1D endif if (size(qsatF ) >= lsizeo) then call mct_aVect_exportRattr(avo,'qsat' ,tempR1D,lsizeo) qsatF = tempR1D endif deallocate(tempR1D) call mct_aVect_clean(avi) call mct_aVect_clean(avo) call mct_gGrid_clean(ggridi) call mct_gsmap_clean(gsmapi) end subroutine datm_shr_getFactors !=============================================================================== real(R8) function datm_shr_eSat(tK,tKbot) !--- arguments --- real(R8),intent(in) :: tK ! temp used in polynomial calculation real(R8),intent(in) :: tKbot ! bottom atm temp !--- local --- real(R8) :: t ! tK converted to Celcius real(R8),parameter :: tkFrz = SHR_CONST_TKFRZ ! freezing T of fresh water ~ K !--- coefficients for esat over water --- real(R8),parameter :: a0=6.107799961_R8 real(R8),parameter :: a1=4.436518521e-01_R8 real(R8),parameter :: a2=1.428945805e-02_R8 real(R8),parameter :: a3=2.650648471e-04_R8 real(R8),parameter :: a4=3.031240396e-06_R8 real(R8),parameter :: a5=2.034080948e-08_R8 real(R8),parameter :: a6=6.136820929e-11_R8 !--- coefficients for esat over ice --- real(R8),parameter :: b0=6.109177956_R8 real(R8),parameter :: b1=5.034698970e-01_R8 real(R8),parameter :: b2=1.886013408e-02_R8 real(R8),parameter :: b3=4.176223716e-04_R8 real(R8),parameter :: b4=5.824720280e-06_R8 real(R8),parameter :: b5=4.838803174e-08_R8 real(R8),parameter :: b6=1.838826904e-10_R8 !---------------------------------------------------------------------------- ! use polynomials to calculate saturation vapor pressure and derivative with ! respect to temperature: over water when t > 0 c and over ice when t <= 0 c ! required to convert relative humidity to specific humidity !---------------------------------------------------------------------------- t = min( 50.0_R8, max(-50.0_R8,(tK-tKfrz)) ) if ( tKbot < tKfrz) then datm_shr_eSat = 100.0_R8*(b0+t*(b1+t*(b2+t*(b3+t*(b4+t*(b5+t*b6)))))) else datm_shr_eSat = 100.0_R8*(a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6)))))) end if end function datm_shr_eSat !=============================================================================== !=============================================================================== end module datm_shr_mod