#ifdef AIX
@PROCESS ALIAS_SIZE(805306368)
#endif
module dice_comp_mod 1,16
! !USES:
use shr_const_mod
use shr_sys_mod
use shr_kind_mod , only: IN=>SHR_KIND_IN, R8=>SHR_KIND_R8, &
CS=>SHR_KIND_CS, CL=>SHR_KIND_CL
use shr_file_mod , only: shr_file_getunit, shr_file_getlogunit, shr_file_getloglevel, &
shr_file_setlogunit, shr_file_setloglevel, shr_file_setio, &
shr_file_freeunit
use shr_mpi_mod , only: shr_mpi_bcast
use shr_flux_mod , only: shr_flux_atmIce
use shr_cal_mod , only: shr_cal_ymd2eday
use mct_mod
use esmf_mod
use perf_mod
use shr_strdata_mod
use shr_dmodel_mod
use shr_pcdf_mod
use seq_cdata_mod
use seq_infodata_mod
use seq_timemgr_mod
use seq_flds_indices , only: nflds_i2x, &
nflds_x2i
use seq_flds_mod , only: seq_flds_i2x_fields, &
seq_flds_x2i_fields
!
! !PUBLIC TYPES:
implicit none
private ! except
!--------------------------------------------------------------------------
! Public interfaces
!--------------------------------------------------------------------------
public :: dice_comp_init
public :: dice_comp_run
public :: dice_comp_final
!--------------------------------------------------------------------------
! Private data
!--------------------------------------------------------------------------
!--- other ---
character(CS) :: myModelName = 'ice' ! user defined model name
integer(IN) :: mpicom
integer(IN) :: my_task ! my task in mpi communicator mpicom
integer(IN) :: npes ! total number of tasks
integer(IN),parameter :: master_task=0 ! task number of master task
integer(IN) :: logunit ! logging unit number
character(CL) :: ice_mode ! mode
integer(IN) :: dbug = 0 ! debug level (higher is more)
logical :: firstcall ! first call logical
logical :: scmMode = .false. ! single column mode
real(R8) :: scmLat = shr_const_SPVAL ! single column lat
real(R8) :: scmLon = shr_const_SPVAL ! single column lon
logical :: read_restart ! start from restart
real(R8) :: flux_swpf ! short-wave penatration factor
real(R8) :: flux_Qmin ! bound on melt rate
logical :: flux_Qacc ! activates water accumulation/melt wrt Q
real(R8) :: flux_Qacc0 ! initial water accumulation value
character(len=*),parameter :: rpfile = 'rpointer.ice'
character(len=*),parameter :: nullstr = 'undefined'
real(R8),parameter :: pi = shr_const_pi ! pi
real(R8),parameter :: spval = shr_const_spval ! flags invalid data
real(R8),parameter :: tFrz = shr_const_tkfrzsw ! temp of freezing salt-water
real(R8),parameter :: latice = shr_const_latice ! latent heat of fusion
real(R8),parameter :: cDay = shr_const_cDay ! sec in calendar day
real(R8),parameter :: waterMax = 1000.0_R8 ! wrt iFrac comp & frazil ice (kg/m^2)
!----- surface albedo constants ------
real(R8),parameter :: snwfrac = 0.286_R8 ! snow cover fraction ~ [0,1]
real(R8),parameter :: as_nidf = 0.950_R8 ! albedo: snow,near-infr,diffuse
real(R8),parameter :: as_vsdf = 0.700_R8 ! albedo: snow,visible ,diffuse
real(R8),parameter :: as_nidr = 0.960_R8 ! albedo: snow,near-infr,direct
real(R8),parameter :: as_vsdr = 0.800_R8 ! albedo: snow,visible ,direct
real(R8),parameter :: ai_nidf = 0.700_R8 ! albedo: ice, near-infr,diffuse
real(R8),parameter :: ai_vsdf = 0.500_R8 ! albedo: ice, visible ,diffuse
real(R8),parameter :: ai_nidr = 0.700_R8 ! albedo: ice, near-infr,direct
real(R8),parameter :: ai_vsdr = 0.500_R8 ! albedo: ice, visible ,direct
real(R8),parameter :: ax_nidf = ai_nidf*(1.0_R8-snwfrac) + as_nidf*snwfrac
real(R8),parameter :: ax_vsdf = ai_vsdf*(1.0_R8-snwfrac) + as_vsdf*snwfrac
real(R8),parameter :: ax_nidr = ai_nidr*(1.0_R8-snwfrac) + as_nidr*snwfrac
real(R8),parameter :: ax_vsdr = ai_vsdr*(1.0_R8-snwfrac) + as_vsdr*snwfrac
integer(IN) :: kswvdr,kswndr,kswvdf,kswndf,kq,kz,kua,kva,kptem,kshum,kdens,ktbot
integer(IN) :: kiFrac,kt,kavsdr,kanidr,kavsdf,kanidf,kswnet,kmelth,kmeltw
integer(IN) :: ksen,klat,klwup,kevap,ktauxa,ktauya,ktref,kqref,kswpen,ktauxo,ktauyo,ksalt
type(shr_strdata_type) :: SDICE
type(mct_rearr) :: rearr
! type(mct_avect) :: avstrm ! av of data from stream
integer(IN) , pointer :: imask(:)
real(R8) , pointer :: yc(:)
real(R8) , pointer :: water(:)
! real(R8) , pointer :: ifrac0(:)
integer(IN),parameter :: ktrans = 42
character(16),parameter :: avofld(1:ktrans) = &
(/"So_t ","So_s ","So_u ","So_v ", &
"So_dhdx ","So_dhdy ","Fioo_q ","Sa_z ", &
"Sa_u ","Sa_v ","Sa_ptem ","Sa_tbot ", &
"Sa_shum ","Sa_dens ","Faxa_swndr ","Faxa_swvdr ", &
"Faxa_swndf ","Faxa_swvdf ","Faxa_lwdn ","Faxa_rain ", &
"Faxa_snow ","Si_t ","Si_tref ","Si_qref ", &
"Si_ifrac ","Si_avsdr ","Si_anidr ","Si_avsdf ", &
"Si_anidf ","Faii_taux ","Faii_tauy ","Faii_lat ", &
"Faii_sen ","Faii_lwup ","Faii_evap ","Faii_swnet ", &
"Fioi_swpen ","Fioi_melth ","Fioi_meltw ","Fioi_salt ", &
"Fioi_taux ","Fioi_tauy " /)
character(16),parameter :: avifld(1:ktrans) = &
(/"to ","s ","uo ","vo ", &
"dhdx ","dhdy ","q ","z ", &
"ua ","va ","ptem ","tbot ", &
"shum ","dens ","swndr ","swvdr ", &
"swndf ","swvdf ","lwdn ","rain ", &
"snow ","t ","tref ","qref ", &
"ifrac ","avsdr ","anidr ","avsdf ", &
"anidf ","tauxa ","tauya ","lat ", &
"sen ","lwup ","evap ","swnet ", &
"swpen ","melth ","meltw ","salt ", &
"tauxo ","tauyo " /)
save
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CONTAINS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: dice_comp_init
!
! !DESCRIPTION:
! initialize data ice model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
subroutine dice_comp_init( EClock, cdata, x2i, i2x, NLFilename ) 2,35
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(ESMF_Clock) , intent(in) :: EClock
type(seq_cdata) , intent(inout) :: cdata
type(mct_aVect) , intent(inout) :: x2i, i2x
character(len=*), optional , intent(in) :: NLFilename ! Namelist filename
!EOP
!--- local variables ---
integer(IN) :: n,k ! generic counters
integer(IN) :: ierr ! error code
integer(IN) :: COMPID ! comp id
integer(IN) :: gsize ! global size
integer(IN) :: lsize ! local size
integer(IN) :: shrlogunit, shrloglev ! original log unit and level
integer(IN) :: nunit ! unit number
integer(IN) :: kfld ! field reference
logical :: ice_present ! flag
logical :: ice_prognostic ! flag
type(seq_infodata_type), pointer :: infodata
type(mct_gsMap) , pointer :: gsmap
type(mct_gGrid) , pointer :: ggrid
character(CL) :: filePath ! generic file path
character(CL) :: fileName ! generic file name
character(CS) :: timeName ! domain file: time variable name
character(CS) :: lonName ! domain file: lon variable name
character(CS) :: latName ! domain file: lat variable name
character(CS) :: maskName ! domain file: mask variable name
character(CS) :: areaName ! domain file: area variable name
integer(IN) :: yearFirst ! first year to use in data stream
integer(IN) :: yearLast ! last year to use in data stream
integer(IN) :: yearAlign ! data year that aligns with yearFirst
character(CL) :: ice_in ! dshr ice namelist
character(CL) :: decomp ! decomp strategy
character(CL) :: rest_file ! restart filename
character(CL) :: rest_file_strm ! restart filename for stream
character(CL) :: restfilm ! model restart file namelist
character(CL) :: restfils ! stream restart file namelist
logical :: exists ! file existance logical
integer(IN) :: nu ! unit number
!----- define namelist -----
namelist / dice_nml / &
ice_in, decomp, flux_swpf, flux_Qmin, flux_Qacc, flux_Qacc0, restfilm, restfils
!--- formats ---
character(*), parameter :: F00 = "('(dice_comp_init) ',8a)"
character(*), parameter :: F01 = "('(dice_comp_init) ',a,5i8)"
character(*), parameter :: F02 = "('(dice_comp_init) ',a,4es13.6)"
character(*), parameter :: F03 = "('(dice_comp_init) ',a,i8,a)"
character(*), parameter :: F04 = "('(dice_comp_init) ',2a,2i8,'s')"
character(*), parameter :: F05 = "('(dice_comp_init) ',a,2f10.4)"
character(*), parameter :: F06 = "('(dice_comp_init) ',a,5l3)"
character(*), parameter :: F90 = "('(dice_comp_init) ',73('='))"
character(*), parameter :: F91 = "('(dice_comp_init) ',73('-'))"
character(*), parameter :: subName = "(dice_comp_init) "
!-------------------------------------------------------------------------------
call t_startf('DICE_INIT')
firstcall = .true.
! Set cdata pointers
call seq_cdata_setptrs(cdata, ID=COMPID, mpicom=mpicom, &
gsMap=gsmap, dom=ggrid, infodata=infodata)
! Determine communicator groups and sizes
call mpi_comm_rank(mpicom, my_task, ierr)
call mpi_comm_size(mpicom, npes, ierr)
!--- open log file ---
if (my_task == master_task) then
logUnit = shr_file_getUnit()
call shr_file_setIO('ice_modelio.nml',logUnit)
else
logUnit = 6
endif
!----------------------------------------------------------------------------
! Reset shr logging to my log file
!----------------------------------------------------------------------------
call shr_file_getLogUnit (shrlogunit)
call shr_file_getLogLevel(shrloglev)
call shr_file_setLogUnit (logUnit)
!----------------------------------------------------------------------------
! Set a Few Defaults
!----------------------------------------------------------------------------
call seq_infodata_getData(infodata,single_column=scmMode, &
& scmlat=scmlat, scmlon=scmLon)
ice_present = .false.
ice_prognostic = .false.
call seq_infodata_GetData(infodata,read_restart=read_restart)
!----------------------------------------------------------------------------
! Read dice_in
!----------------------------------------------------------------------------
call t_startf('dice_readnml')
filename = "dice_in"
ice_in = "unset"
decomp = "1d"
flux_swpf = 0.0_R8 ! no penetration
flux_Qmin = -300.0_R8 ! kg/s/m^2
flux_Qacc = .false. ! no accumulation
flux_Qacc0 = 0.0_R8 ! no water
restfilm = trim(nullstr)
restfils = trim(nullstr)
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=dice_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)' ice_in = ',trim(ice_in)
write(logunit,F00)' decomp = ',trim(decomp)
write(logunit,F02)' flux_swpf = ',flux_swpf
write(logunit,F02)' flux_Qmin = ',flux_Qmin
write(logunit,F06)' flux_Qacc = ',flux_Qacc
write(logunit,F02)' flux_Qacc0 = ',flux_Qacc0
write(logunit,F00)' restfilm = ',trim(restfilm)
write(logunit,F00)' restfils = ',trim(restfils)
endif
call shr_mpi_bcast(ice_in ,mpicom,'ice_in')
call shr_mpi_bcast(decomp ,mpicom,'decomp')
call shr_mpi_bcast(flux_swpf ,mpicom,'flux_swpf')
call shr_mpi_bcast(flux_Qmin ,mpicom,'flux_Qmin')
call shr_mpi_bcast(flux_Qacc ,mpicom,'flux_Qacc')
call shr_mpi_bcast(flux_Qacc0,mpicom,'flux_Qacc0')
call shr_mpi_bcast(restfilm,mpicom,'restfilm')
call shr_mpi_bcast(restfils,mpicom,'restfils')
rest_file = trim(restfilm)
rest_file_strm = trim(restfils)
!----------------------------------------------------------------------------
! Read dshr namelist
!----------------------------------------------------------------------------
call shr_strdata_readnml(SDICE,trim(ice_in),mpicom=mpicom)
!----------------------------------------------------------------------------
! Validate mode
!----------------------------------------------------------------------------
ice_mode = trim(SDICE%dataMode)
! check that we know how to handle the mode
if (trim(ice_mode) == 'NULL' .or. &
trim(ice_mode) == 'SSTDATA' .or. &
trim(ice_mode) == 'COPYALL') then
if (my_task == master_task) &
write(logunit,F00) ' ice mode = ',trim(ice_mode)
else
write(logunit,F00) ' ERROR illegal ice mode = ',trim(ice_mode)
call shr_sys_abort()
endif
call t_stopf('dice_readnml')
!----------------------------------------------------------------------------
! Initialize datasets
!----------------------------------------------------------------------------
call t_startf('dice_strdata_init')
if (trim(ice_mode) /= 'NULL') then
ice_present = .true.
if (scmmode) then
if (my_task == master_task) &
write(logunit,F05) ' scm lon lat = ',scmlon,scmlat
call shr_strdata_init(SDICE,mpicom,compid,name='ice', &
scmmode=scmmode,scmlon=scmlon,scmlat=scmlat)
else
call shr_strdata_init(SDICE,mpicom,compid,name='ice')
endif
endif
if (trim(ice_mode) == 'SSTDATA' .or. &
trim(ice_mode) == 'COPYALL') then
ice_prognostic = .true.
endif
if (my_task == master_task) then
call shr_strdata_print(SDICE,'SDICE data')
endif
call t_stopf('dice_strdata_init')
!----------------------------------------------------------------------------
! Set flag to specify data components
!----------------------------------------------------------------------------
call seq_infodata_PutData(infodata, &
ice_present=ice_present, ice_prognostic=ice_prognostic, &
ice_nx=SDICE%nxg, ice_ny=SDICE%nyg )
!----------------------------------------------------------------------------
! Initialize MCT global seg map, 1d decomp
!----------------------------------------------------------------------------
call t_startf('dice_initgsmaps')
if (my_task == master_task) write(logunit,F00) ' initialize gsmaps'
call shr_sys_flush(logunit)
call shr_dmodel_gsmapcreate(gsmap,SDICE%nxg*SDICE%nyg,compid,mpicom,decomp)
lsize = mct_gsmap_lsize(gsmap,mpicom)
if (ice_present) then
call mct_rearr_init(SDICE%gsmap,gsmap,mpicom,rearr)
endif
call t_stopf('dice_initgsmaps')
!----------------------------------------------------------------------------
! Initialize MCT domain
!----------------------------------------------------------------------------
call t_startf('dice_initmctdom')
if (my_task == master_task) write(logunit,F00) 'copy domains'
call shr_sys_flush(logunit)
if (ice_present) call shr_dmodel_rearrGGrid(SDICE%grid, ggrid, gsmap, rearr, mpicom)
call t_stopf('dice_initmctdom')
!----------------------------------------------------------------------------
! Initialize MCT attribute vectors
!----------------------------------------------------------------------------
call t_startf('dice_initmctavs')
if (my_task == master_task) write(logunit,F00) 'allocate AVs'
call shr_sys_flush(logunit)
call mct_aVect_init(i2x, rList=seq_flds_i2x_fields, lsize=lsize)
call mct_aVect_zero(i2x)
kiFrac = mct_aVect_indexRA(i2x,'Si_ifrac')
kt = mct_aVect_indexRA(i2x,'Si_t')
ktref = mct_aVect_indexRA(i2x,'Si_tref')
kqref = mct_aVect_indexRA(i2x,'Si_qref')
kavsdr = mct_aVect_indexRA(i2x,'Si_avsdr')
kanidr = mct_aVect_indexRA(i2x,'Si_anidr')
kavsdf = mct_aVect_indexRA(i2x,'Si_avsdf')
kanidf = mct_aVect_indexRA(i2x,'Si_anidf')
kswnet = mct_aVect_indexRA(i2x,'Faii_swnet')
ksen = mct_aVect_indexRA(i2x,'Faii_sen')
klat = mct_aVect_indexRA(i2x,'Faii_lat')
klwup = mct_aVect_indexRA(i2x,'Faii_lwup')
kevap = mct_aVect_indexRA(i2x,'Faii_evap')
ktauxa = mct_aVect_indexRA(i2x,'Faii_taux')
ktauya = mct_aVect_indexRA(i2x,'Faii_tauy')
kmelth = mct_aVect_indexRA(i2x,'Fioi_melth')
kmeltw = mct_aVect_indexRA(i2x,'Fioi_meltw')
kswpen = mct_aVect_indexRA(i2x,'Fioi_swpen')
ktauxo = mct_aVect_indexRA(i2x,'Fioi_taux')
ktauyo = mct_aVect_indexRA(i2x,'Fioi_tauy')
ksalt = mct_aVect_indexRA(i2x,'Fioi_salt')
call mct_aVect_init(x2i, rList=seq_flds_x2i_fields, lsize=lsize)
call mct_aVect_zero(x2i)
kswvdr = mct_aVect_indexRA(x2i,'Faxa_swvdr')
kswndr = mct_aVect_indexRA(x2i,'Faxa_swndr')
kswvdf = mct_aVect_indexRA(x2i,'Faxa_swvdf')
kswndf = mct_aVect_indexRA(x2i,'Faxa_swndf')
kq = mct_aVect_indexRA(x2i,'Fioo_q')
kz = mct_aVect_indexRA(x2i,'Sa_z')
kua = mct_aVect_indexRA(x2i,'Sa_u')
kva = mct_aVect_indexRA(x2i,'Sa_v')
kptem = mct_aVect_indexRA(x2i,'Sa_ptem')
kshum = mct_aVect_indexRA(x2i,'Sa_shum')
kdens = mct_aVect_indexRA(x2i,'Sa_dens')
ktbot = mct_aVect_indexRA(x2i,'Sa_tbot')
! call mct_aVect_init(avstrm, rList=flds_strm, lsize=lsize)
! call mct_aVect_zero(avstrm)
allocate(imask(lsize))
allocate(yc(lsize))
allocate(water(lsize))
! allocate(iFrac0(lsize))
kfld = mct_aVect_indexRA(ggrid%data,'mask')
imask(:) = nint(ggrid%data%rAttr(kfld,:))
kfld = mct_aVect_indexRA(ggrid%data,'lat')
yc(:) = ggrid%data%rAttr(kfld,:)
call t_stopf('dice_initmctavs')
!----------------------------------------------------------------------------
! Read restart
!----------------------------------------------------------------------------
if (read_restart) then
if (trim(rest_file) == trim(nullstr) .and. &
trim(rest_file_strm) == trim(nullstr)) then
if (my_task == master_task) then
write(logunit,F00) ' restart filenames from rpointer'
call shr_sys_flush(logunit)
inquire(file=trim(rpfile),exist=exists)
if (.not.exists) then
write(logunit,F00) ' ERROR: rpointer file does not exist'
call shr_sys_abort(trim(subname)//' ERROR: rpointer file missing')
endif
nu = shr_file_getUnit()
open(nu,file=trim(rpfile),form='formatted')
read(nu,'(a)') rest_file
read(nu,'(a)') rest_file_strm
close(nu)
call shr_file_freeUnit(nu)
inquire(file=trim(rest_file_strm),exist=exists)
endif
call shr_mpi_bcast(rest_file,mpicom,'rest_file')
call shr_mpi_bcast(rest_file_strm,mpicom,'rest_file_strm')
else
! use namelist already read
if (my_task == master_task) then
write(logunit,F00) ' restart filenames from namelist '
call shr_sys_flush(logunit)
inquire(file=trim(rest_file_strm),exist=exists)
endif
endif
call shr_mpi_bcast(exists,mpicom,'exists')
if (my_task == master_task) write(logunit,F00) ' reading ',trim(rest_file)
call shr_pcdf_readwrite('read',trim(rest_file),mpicom,gsmap,rf1=water,rf1n='water')
if (exists) then
if (my_task == master_task) write(logunit,F00) ' reading ',trim(rest_file_strm)
call shr_strdata_restRead(trim(rest_file_strm),SDICE,mpicom)
else
if (my_task == master_task) write(logunit,F00) ' file not found, skipping ',trim(rest_file_strm)
endif
call shr_sys_flush(logunit)
endif
!----------------------------------------------------------------------------
! On initial call, x2i is unset, so set for use in run method
! These values should have no impact on the solution!!
!----------------------------------------------------------------------------
x2i%rAttr(kz,:) = 10.0_R8
x2i%rAttr(kua,:) = 5.0_R8
x2i%rAttr(kva,:) = 5.0_R8
x2i%rAttr(kptem,:) = 260.0_R8
x2i%rAttr(ktbot,:) = 260.0_R8
x2i%rAttr(kshum,:) = 0.0014_R8
x2i%rAttr(kdens,:) = 1.3_R8
!----------------------------------------------------------------------------
! Set initial ice state, needed for CCSM atm initialization
!----------------------------------------------------------------------------
call t_adj_detailf(+2)
call dice_comp_run( EClock, cdata, x2i, i2x)
call t_adj_detailf(-2)
!----------------------------------------------------------------------------
! Reset shr logging to original values
!----------------------------------------------------------------------------
if (my_task == master_task) write(logunit,F00) 'dice_comp_init done'
call shr_sys_flush(logunit)
call shr_file_setLogUnit (shrlogunit)
call shr_file_setLogLevel(shrloglev)
call shr_sys_flush(logunit)
call t_stopf('DICE_INIT')
end subroutine dice_comp_init
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: dice_comp_run
!
! !DESCRIPTION:
! run method for dead ice model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
subroutine dice_comp_run( EClock, cdata, x2i, i2x) 2,24
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(ESMF_Clock) ,intent(in) :: EClock
type(seq_cdata) ,intent(inout) :: cdata
type(mct_aVect) ,intent(inout) :: x2i ! driver -> dead
type(mct_aVect) ,intent(inout) :: i2x ! dead -> driver
!EOP
!--- local ---
type(mct_gsMap) , pointer :: gsmap
type(mct_gGrid) , pointer :: ggrid
integer(IN) :: CurrentYMD ! model date
integer(IN) :: CurrentTOD ! model sec into model date
integer(IN) :: yy,mm,dd ! year month day
integer(IN) :: n ! indices
integer(IN) :: nf ! fields loop index
integer(IN) :: nl ! ice frac index
integer(IN) :: lsize ! size of attr vect
integer(IN) :: shrlogunit, shrloglev ! original log unit and level
logical :: glcrun_alarm ! is glc going to run now
logical :: newdata ! has newdata been read
logical :: mssrmlf ! remove old data
integer(IN) :: idt ! integer timestep
real(R8) :: dt ! timestep
real(R8) :: hn ! h field
logical :: write_restart ! restart now
character(CL) :: case_name ! case name
character(CL) :: rest_file ! restart_file
character(CL) :: rest_file_strm ! restart_file for stream
integer(IN) :: nu ! unit number
real(R8) :: qmeltall ! q that would melt all accumulated water
real(R8) :: cosarg ! for setting ice temp pattern
integer(IN) :: eday, eday0 ! elapsed day counters
type(seq_infodata_type), pointer :: infodata
character(*), parameter :: F00 = "('(dice_comp_run) ',8a)"
character(*), parameter :: F04 = "('(dice_comp_run) ',2a,2i8,'s')"
character(*), parameter :: subName = "(dice_comp_run) "
!-------------------------------------------------------------------------------
call t_startf('DICE_RUN')
call t_startf('dice_run1')
!----------------------------------------------------------------------------
! Reset shr logging to my log file
!----------------------------------------------------------------------------
call shr_file_getLogUnit (shrlogunit)
call shr_file_getLogLevel(shrloglev)
call shr_file_setLogUnit (logUnit)
call seq_cdata_setptrs(cdata, gsMap=gsmap, dom=ggrid)
call seq_cdata_setptrs(cdata, infodata=infodata)
call seq_timemgr_EClockGetData( EClock, curr_ymd=CurrentYMD, curr_tod=CurrentTOD)
call seq_timemgr_EClockGetData( EClock, curr_yr=yy, curr_mon=mm, curr_day=dd)
call seq_timemgr_EClockGetData( EClock, dtime=idt)
dt = idt * 1.0_r8
write_restart = seq_timemgr_RestartAlarmIsOn(EClock)
call t_stopf('dice_run1')
!--------------------
! UNPACK
!--------------------
call t_startf('dice_unpack')
! lsize = mct_avect_lsize(x2i)
! do nf=1,nflds_x2i
! do n=1,lsize
! ?? = x2i%rAttr(nf,n)
! enddo
! enddo
call t_stopf('dice_unpack')
!--------------------
! ADVANCE ICE
!--------------------
call t_barrierf('dice_BARRIER',mpicom)
call t_startf('dice')
!--- copy all fields from streams to i2x as default ---
if (trim(ice_mode) /= 'NULL') then
call t_startf('dice_strdata_advance')
call shr_strdata_advance(SDICE,currentYMD,currentTOD,mpicom,'dice')
call t_stopf('dice_strdata_advance')
call t_barrierf('dice_scatter_BARRIER',mpicom)
call t_startf('dice_scatter')
do n = 1,SDICE%nstreams
call shr_dmodel_translateAV(SDICE%avs(n),i2x,avifld,avofld,rearr)
enddo
call t_stopf('dice_scatter')
else
call mct_aVect_zero(i2x)
endif
call t_startf('dice_mode')
select case (trim(ice_mode))
case('COPYALL')
! do nothing extra
case('SSTDATA')
if (firstcall .and. .not. read_restart) then
! iFrac0 = iFrac ! previous step's ice fraction
water = 0.0_R8 ! previous step's water accumulation
where (i2x%rAttr(kiFrac,:) > 0.0_R8) water(:) = flux_Qacc0
endif
call shr_cal_ymd2eday(0,mm,dd,eDay ) ! model date
call shr_cal_ymd2eday(0,09,01,eDay0) ! sept 1st
cosArg = 2.0_R8*pi*(real(eDay,R8) + real(currentTOD,R8)/cDay - real(eDay0,R8))/365.0_R8
lsize = mct_avect_lsize(i2x)
do n = 1,lsize
!--- fix erroneous iFrac ---
i2x%rAttr(kiFrac,n) = min(1.0_R8,max(0.0_R8,i2x%rAttr(kiFrac,n)))
!--- fabricate ice surface T, fix erroneous iFrac ---
if ( yc(n) > 0.0_R8) then
i2x%rAttr(kt,n) = 260.0_R8 + 10.0_R8*cos(cosArg)
else
i2x%rAttr(kt,n) = 260.0_R8 - 10.0_R8*cos(cosArg)
end if
!--- set albedos (constant) ---
i2x%rAttr(kavsdr,n) = ax_vsdr
i2x%rAttr(kanidr,n) = ax_nidr
i2x%rAttr(kavsdf,n) = ax_vsdf
i2x%rAttr(kanidf,n) = ax_nidf
!--- swnet is sent to cpl as a diagnostic quantity only ---
!--- newly recv'd swdn goes with previously sent albedo ---
!--- but albedos are (currently) time invariant ---
i2x%rAttr(kswnet,n) = (1.0_R8 - i2x%rAttr(kavsdr,n))*x2i%rAttr(kswvdr,n) &
& + (1.0_R8 - i2x%rAttr(kanidr,n))*x2i%rAttr(kswndr,n) &
& + (1.0_R8 - i2x%rAttr(kavsdf,n))*x2i%rAttr(kswvdf,n) &
& + (1.0_R8 - i2x%rAttr(kanidf,n))*x2i%rAttr(kswndf,n)
!--- compute melt/freeze water balance, adjust iFrac -------------
if ( .not. flux_Qacc ) then ! Q accumulation option is OFF
i2x%rAttr(kmelth,n) = min(x2i%rAttr(kq,n),0.0_R8 ) ! q<0 => melt potential
i2x%rAttr(kmelth,n) = max(i2x%rAttr(kmelth,n),Flux_Qmin ) ! limit the melt rate
i2x%rAttr(kmeltw,n) = -i2x%rAttr(kmelth,n)/latice ! corresponding water flux
else ! Q accumulation option is ON
!--------------------------------------------------------------
! 1a) Q<0 & iFrac > 0 => infinite supply of water to melt
! 1b) Q<0 & iFrac = 0 => melt accumulated water only
! 2a) Q>0 & iFrac > 0 => zero-out accumulated water
! 2b) Q>0 & iFrac = 0 => accumulated water
!--------------------------------------------------------------
if ( x2i%rAttr(kq,n) < 0.0_R8 ) then ! Q<0 => melt
if (i2x%rAttr(kiFrac,n) > 0.0_R8 ) then
i2x%rAttr(kmelth,n) = i2x%rAttr(kiFrac,n)*max(x2i%rAttr(kq,n),Flux_Qmin)
i2x%rAttr(kmeltw,n) = -i2x%rAttr(kmelth,n)/latice
! water(n) = < don't change this value >
else
Qmeltall = -water(n)*latice/dt
i2x%rAttr(kmelth,n) = max(x2i%rAttr(kq,n), Qmeltall, Flux_Qmin )
i2x%rAttr(kmeltw,n) = -i2x%rAttr(kmelth,n)/latice
water(n) = water(n) - i2x%rAttr(kmeltw,n)*dt
end if
else ! Q>0 => freeze
if (i2x%rAttr(kiFrac,n) > 0.0_R8 ) then
i2x%rAttr(kmelth,n) = 0.0_R8
i2x%rAttr(kmeltw,n) = 0.0_R8
water(n) = 0.0_R8
else
i2x%rAttr(kmelth,n) = 0.0_R8
i2x%rAttr(kmeltw,n) = 0.0_R8
water(n) = water(n) + dt*i2x%rAttr(kq,n)/latice
end if
end if
if (water(n) < 1.0e-16_R8 ) water(n) = 0.0_R8
!--- non-zero water => non-zero iFrac ---
if (i2x%rAttr(kiFrac,n) <= 0.0_R8 .and. water(n) > 0.0_R8) then
i2x%rAttr(kiFrac,n) = min(1.0_R8,water(n)/waterMax)
! i2x%rAttr(kT,n) = Tfrz ! T can be above freezing?!?
end if
!--- cpl multiplies melth & meltw by iFrac ---
!--- divide by iFrac here => fixed quantity flux (not per area) ---
if (i2x%rAttr(kiFrac,n) > 0.0_R8) then
i2x%rAttr(kiFrac,n) = max( 0.01_R8, i2x%rAttr(kiFrac,n)) ! min iFrac
i2x%rAttr(kmelth,n) = i2x%rAttr(kmelth,n)/i2x%rAttr(kiFrac,n)
i2x%rAttr(kmeltw,n) = i2x%rAttr(kmeltw,n)/i2x%rAttr(kiFrac,n)
else
i2x%rAttr(kmelth,n) = 0.0_R8
i2x%rAttr(kmeltw,n) = 0.0_R8
end if
end if
!--- modify T wrt iFrac: (iFrac -> 0) => (T -> Tfrz) ---
i2x%rAttr(kt,n) = Tfrz + i2x%rAttr(kiFrac,n)*(i2x%rAttr(kt,n)-Tfrz)
end do
!----------------------------------------------------------------------------
! compute atm/ice surface fluxes
!----------------------------------------------------------------------------
call shr_flux_atmIce(iMask ,x2i%rAttr(kz,:) ,x2i%rAttr(kua,:) ,x2i%rAttr(kva,:), &
x2i%rAttr(kptem,:) ,x2i%rAttr(kshum,:) ,x2i%rAttr(kdens,:) ,x2i%rAttr(ktbot,:), &
i2x%rAttr(kt,:) ,i2x%rAttr(ksen,:) ,i2x%rAttr(klat,:) ,i2x%rAttr(klwup,:), &
i2x%rAttr(kevap,:) ,i2x%rAttr(ktauxa,:) ,i2x%rAttr(ktauya,:) ,i2x%rAttr(ktref,:), &
i2x%rAttr(kqref,:) )
!----------------------------------------------------------------------------
! compute ice/oce surface fluxes (except melth & meltw, see above)
!----------------------------------------------------------------------------
do n=1,lsize
if (iMask(n) == 0) then
i2x%rAttr(kswpen,n) = spval
i2x%rAttr(kmelth,n) = spval
i2x%rAttr(kmeltw,n) = spval
i2x%rAttr(ksalt ,n) = spval
i2x%rAttr(ktauxo,n) = spval
i2x%rAttr(ktauyo,n) = spval
i2x%rAttr(kiFrac,n) = 0.0_R8
else
!--- penetrating short wave ---
i2x%rAttr(kswpen,n) = max(0.0_R8, flux_swpf*i2x%rAttr(kswnet,n) ) ! must be non-negative
!--- i/o surface stress ( = atm/ice stress) ---
i2x%rAttr(ktauxo,n) = i2x%rAttr(ktauxa,n)
i2x%rAttr(ktauyo,n) = i2x%rAttr(ktauya,n)
!--- salt flux ---
i2x%rAttr(ksalt ,n) = 0.0_R8
end if
! !--- save ifrac for next timestep
! iFrac0(n) = i2x%rAttr(kiFrac,n)
end do
end select
call t_stopf('dice_mode')
if (write_restart) then
call t_startf('dice_restart')
call seq_infodata_GetData( infodata, case_name=case_name)
write(rest_file,"(2a,i4.4,a,i2.2,a,i2.2,a,i5.5,a)") &
trim(case_name), '.dice.r.', yy,'-',mm,'-',dd,'-',currentTOD,'.nc'
write(rest_file_strm,"(2a,i4.4,a,i2.2,a,i2.2,a,i5.5,a)") &
trim(case_name), '.dice.rs1.', yy,'-',mm,'-',dd,'-',currentTOD,'.bin'
if (my_task == master_task) then
nu = shr_file_getUnit()
open(nu,file=trim(rpfile),form='formatted')
write(nu,'(a)') rest_file
write(nu,'(a)') rest_file_strm
close(nu)
call shr_file_freeUnit(nu)
endif
if (my_task == master_task) write(logunit,F04) ' writing ',trim(rest_file),currentYMD,currentTOD
call shr_pcdf_readwrite('write',trim(rest_file),mpicom,gsmap,clobber=.true., &
rf1=water,rf1n='water')
if (my_task == master_task) write(logunit,F04) ' writing ',trim(rest_file_strm),currentYMD,currentTOD
call shr_strdata_restWrite(trim(rest_file_strm),SDICE,mpicom,trim(case_name),'SDICE strdata')
call shr_sys_flush(logunit)
call t_stopf('dice_restart')
endif
call t_stopf('dice')
!----------------------------------------------------------------------------
! Log output for model date
! Reset shr logging to original values
!----------------------------------------------------------------------------
call t_startf('dice_run2')
if (my_task == master_task) then
write(logunit,F04) trim(myModelName),': model date ', CurrentYMD,CurrentTOD
call shr_sys_flush(logunit)
end if
firstcall = .false.
call shr_file_setLogUnit (shrlogunit)
call shr_file_setLogLevel(shrloglev)
call shr_sys_flush(logunit)
call t_stopf('dice_run2')
call t_stopf('DICE_RUN')
end subroutine dice_comp_run
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: dice_comp_final
!
! !DESCRIPTION:
! finalize method for dead ice model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
!
subroutine dice_comp_final() 1
implicit none
!EOP
!--- formats ---
character(*), parameter :: F00 = "('(dice_comp_final) ',8a)"
character(*), parameter :: F91 = "('(dice_comp_final) ',73('-'))"
character(*), parameter :: subName = "(dice_comp_final) "
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
call t_startf('DICE_FINAL')
if (my_task == master_task) then
write(logunit,F91)
write(logunit,F00) trim(myModelName),': end of main integration loop'
write(logunit,F91)
end if
call t_stopf('DICE_FINAL')
end subroutine dice_comp_final
!===============================================================================
!===============================================================================
end module dice_comp_mod