!======================================================================= ! ! Contains Icepack component driver routines common to all drivers. ! ! Authors: Lorenzo Zampieri ( lorenzo.zampieri@awi.de ) ! Modified by Qian Wang to apply to SCHISM !======================================================================= submodule (icedrv_main) icedrv_step !use icedrv_constants, only: c0, c4, nu_diag, ice_stderr use icedrv_kinds use icedrv_system, only: icedrv_system_abort use icepack_intfc, only: icepack_warnings_flush use icepack_intfc, only: icepack_warnings_aborted use icepack_intfc, only: icepack_query_tracer_flags use icepack_intfc, only: icepack_query_tracer_indices use icepack_intfc, only: icepack_query_tracer_sizes use icepack_intfc, only: icepack_query_parameters implicit none !======================================================================= contains !======================================================================= ! ! Scales radiation fields computed on the previous time step. ! subroutine prep_radiation () ! column package includes use icepack_intfc, only: icepack_prep_radiation implicit none ! local variables integer (kind=int_kind) :: & i ! horizontal indices character(len=*), parameter :: subname='(prep_radiation)' !----------------------------------------------------------------- ! Compute netsw scaling factor (new netsw / old netsw) !----------------------------------------------------------------- do i = 1, nx alvdr_init(i) = alvdr_ai(i) alvdf_init(i) = alvdf_ai(i) alidr_init(i) = alidr_ai(i) alidf_init(i) = alidf_ai(i) call icepack_prep_radiation(ncat=ncat, nilyr=nilyr, nslyr=nslyr, & aice=aice(i), aicen=aicen(i,:), & swvdr=swvdr(i), swvdf=swvdf(i), & swidr=swidr(i), swidf=swidf(i), & alvdr_ai=alvdr_ai(i), alvdf_ai=alvdf_ai(i), & alidr_ai=alidr_ai(i), alidf_ai=alidf_ai(i), & scale_factor=scale_factor(i), & fswsfcn=fswsfcn(i,:), fswintn=fswintn(i,:), & fswthrun=fswthrun(i,:), fswpenln=fswpenln(i,:,:), & Sswabsn=Sswabsn(i,:,:), Iswabsn=Iswabsn(i,:,:)) enddo ! i call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) end subroutine prep_radiation !======================================================================= ! ! Driver for updating ice and snow internal temperatures and ! computing thermodynamic growth rates and coupler fluxes. ! subroutine step_therm1 (dt) ! column packge includes use icepack_intfc, only: icepack_step_therm1 use schism_glbl, only : idry use mice_module, only : stress_atmice_x,stress_atmice_y implicit none logical (kind=log_kind) :: & prescribed_ice ! if .true., use prescribed ice instead of computed real (kind=dbl_kind), intent(in) :: & dt ! time step ! local variables integer (kind=int_kind) :: & i , & ! horizontal indices n , & ! thickness category index k, kk ! indices for aerosols integer (kind=int_kind) :: & ntrcr, nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, nt_vlvl, nt_Tsfc, & nt_iage, nt_FY, nt_qice, nt_sice, nt_aero, nt_qsno logical (kind=log_kind) :: & tr_iage, tr_FY, tr_aero, tr_pond, tr_pond_cesm, & tr_pond_lvl, tr_pond_topo, calc_Tsfc, calc_strair real (kind=dbl_kind), dimension(n_aero,2,ncat) :: & aerosno, aeroice ! kg/m^2 real (kind=dbl_kind) :: & puny character(len=*), parameter :: subname='(step_therm1)' !----------------------------------------------------------------- ! query icepack values !----------------------------------------------------------------- call icepack_query_parameters(puny_out=puny) call icepack_query_parameters(calc_strair_out=calc_strair,calc_Tsfc_out=calc_Tsfc) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_tracer_sizes( & ntrcr_out=ntrcr) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_tracer_flags( & tr_iage_out=tr_iage, tr_FY_out=tr_FY, & tr_aero_out=tr_aero, tr_pond_out=tr_pond, tr_pond_cesm_out=tr_pond_cesm, & tr_pond_lvl_out=tr_pond_lvl, tr_pond_topo_out=tr_pond_topo) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_tracer_indices( & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, nt_Tsfc_out=nt_Tsfc, & nt_iage_out=nt_iage, nt_FY_out=nt_FY, & nt_qice_out=nt_qice, nt_sice_out=nt_sice, & nt_aero_out=nt_aero, nt_qsno_out=nt_qsno) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) !----------------------------------------------------------------- prescribed_ice = .false. aerosno(:,:,:) = c0 aeroice(:,:,:) = c0 do i = 1, nx !if(idry(i)==1) cycle !----------------------------------------------------------------- ! Save the ice area passed to the coupler (so that history fields ! can be made consistent with coupler fields). ! Save the initial ice area and volume in each category. !----------------------------------------------------------------- aice_init (i) = aice (i) do n = 1, ncat aicen_init(i,n) = aicen(i,n) vicen_init(i,n) = vicen(i,n) vsnon_init(i,n) = vsnon(i,n) enddo enddo ! i do i = 1, nx !if(idry(i)==1) cycle if (tr_aero) then ! trcrn(nt_aero) has units kg/m^3 do n=1,ncat do k=1,n_aero aerosno (k,:,n) = & trcrn(i,nt_aero+(k-1)*4 :nt_aero+(k-1)*4+1,n) & * vsnon_init(i,n) aeroice (k,:,n) = & trcrn(i,nt_aero+(k-1)*4+2:nt_aero+(k-1)*4+3,n) & * vicen_init(i,n) enddo enddo endif ! tr_aero call icepack_step_therm1(dt=dt, ncat=ncat, nilyr=nilyr, nslyr=nslyr, & aicen_init = aicen_init(i,:), & vicen_init = vicen_init(i,:), & vsnon_init = vsnon_init(i,:), & aice = aice(i), aicen = aicen(i,:), & vice = vice(i), vicen = vicen(i,:), & vsno = vsno(i), vsnon = vsnon(i,:), & uvel = uvel(i), vvel = vvel(i), & Tsfc = trcrn(i,nt_Tsfc,:), & zqsn = trcrn(i,nt_qsno:nt_qsno+nslyr-1,:), & zqin = trcrn(i,nt_qice:nt_qice+nilyr-1,:), & zSin = trcrn(i,nt_sice:nt_sice+nilyr-1,:), & alvl = trcrn(i,nt_alvl,:), & vlvl = trcrn(i,nt_vlvl,:), & apnd = trcrn(i,nt_apnd,:), & hpnd = trcrn(i,nt_hpnd,:), & ipnd = trcrn(i,nt_ipnd,:), & iage = trcrn(i,nt_iage,:), & FY = trcrn(i,nt_FY,:), & aerosno = aerosno(:,:,:), & aeroice = aeroice(:,:,:), & uatm = uatm(i), vatm = vatm(i), & wind = wind(i), & zlvl = zlvl_v, & zlvs = zlvs, & !zlvl_t = zlvl_t, & Qa = Qa(i), rhoa = rhoa(i), & Tair = T_air(i), Tref = Tref(i), & Qref = Qref(i), Uref = Uref(i), & Cdn_atm_ratio = Cdn_atm_ratio(i),& Cdn_ocn = Cdn_ocn(i), & Cdn_ocn_skin = Cdn_ocn_skin(i), & Cdn_ocn_floe = Cdn_ocn_floe(i), & Cdn_ocn_keel = Cdn_ocn_keel(i), & Cdn_atm = Cdn_atm(i), & Cdn_atm_skin = Cdn_atm_skin(i), & Cdn_atm_floe = Cdn_atm_floe(i), & Cdn_atm_pond = Cdn_atm_pond(i), & Cdn_atm_rdg = Cdn_atm_rdg(i), & hfreebd = hfreebd(i), hkeel = hkeel(i), & hdraft = hdraft(i), hridge = hridge(i), & distrdg = distrdg(i), dkeel = dkeel(i), & lfloe = lfloe(i), dfloe = dfloe(i), & strax = strax(i), stray = stray(i), & strairxT = strairxT(i), strairyT = strairyT(i), & potT = potT(i), sst = sst(i), & sss = sss(i), Tf = Tf(i), & strocnxT = strocnxT(i), strocnyT = strocnyT(i), & fbot = fbot(i), frzmlt = frzmlt(i), & Tbot = Tbot(i), Tsnice = Tsnice(i), & rside = rside(i), fside = fside(i), & fsnow = fsnow(i), frain = frain(i), & fpond = fpond(i), & fsurf = fsurf(i), fsurfn = fsurfn(i,:), & fcondtop = fcondtop(i), fcondtopn = fcondtopn(i,:), & fcondbot = fcondbot(i), fcondbotn = fcondbotn(i,:), & fswsfcn = fswsfcn(i,:), fswintn = fswintn(i,:), & fswthrun = fswthrun(i,:), fswabs = fswabs(i), & flwout = flwout(i), flw = flw(i), & fsens = fsens(i), fsensn = fsensn(i,:), & flat = flat(i), flatn = flatn(i,:), & fresh = fresh(i), fsalt = fsalt(i), & fhocn = fhocn(i), fswthru = fswthru(i), & flatn_f = flatn_f(i,:), fsensn_f = fsensn_f(i,:), & fsurfn_f = fsurfn_f(i,:), & fcondtopn_f = fcondtopn_f(i,:), & faero_atm = faero_atm(i,1:n_aero), & faero_ocn = faero_ocn(i,1:n_aero), & Sswabsn = Sswabsn(i,:,:),Iswabsn = Iswabsn(i,:,:), & evap = evap(i), evaps = evaps(i), evapi = evapi(i), & dhsn = dhsn(i,:), ffracn = ffracn(i,:), & meltt = meltt(i), melttn = melttn(i,:), & meltb = meltb(i), meltbn = meltbn(i,:), & melts = melts(i), meltsn = meltsn(i,:), & congel = congel(i), congeln = congeln(i,:), & snoice = snoice(i), snoicen = snoicen(i,:), & dsnown = dsnown(i,:), & lmask_n = lmask_n(i), lmask_s = lmask_s(i), & mlt_onset=mlt_onset(i), frz_onset = frz_onset(i), & yday = yday, prescribed_ice = prescribed_ice) if (calc_strair) then stress_atmice_x(i) = strairxT(i) stress_atmice_y(i) = strairyT(i) endif if (tr_aero) then do n = 1, ncat if (vicen(i,n) > puny) & aeroice(:,:,n) = aeroice(:,:,n)/vicen(i,n) if (vsnon(i,n) > puny) & aerosno(:,:,n) = aerosno(:,:,n)/vsnon(i,n) do k = 1, n_aero do kk = 1, 2 trcrn(i,nt_aero+(k-1)*4+kk-1,n)=aerosno(k,kk,n) trcrn(i,nt_aero+(k-1)*4+kk+1,n)=aeroice(k,kk,n) enddo enddo enddo endif ! tr_aero enddo ! i call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) end subroutine step_therm1 !======================================================================= ! Driver for thermodynamic changes not needed for coupling: ! transport in thickness space, lateral growth and melting. ! subroutine step_therm2 (dt) ! column package_includes use icepack_intfc, only: icepack_step_therm2 use schism_glbl, only : idry implicit none real (kind=dbl_kind), intent(in) :: & dt ! time step ! local variables integer (kind=int_kind) :: & i ! horizontal index integer (kind=int_kind) :: & ntrcr, nbtrcr logical (kind=log_kind) :: & tr_fsd ! floe size distribution tracers character(len=*), parameter :: subname='(step_therm2)' !----------------------------------------------------------------- ! query icepack values !----------------------------------------------------------------- call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_out=nbtrcr) call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) !----------------------------------------------------------------- do i = 1, nx !if(idry(i)==1) cycle ! wave_sig_ht - compute here to pass to add new ice if (tr_fsd) & wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(:))) call icepack_step_therm2(dt=dt, ncat=ncat, & nltrcr=nltrcr, nilyr=nilyr, nslyr=nslyr, & hin_max=hin_max(:), nblyr=nblyr, & aicen=aicen(i,:), & vicen=vicen(i,:), & vsnon=vsnon(i,:), & aicen_init=aicen_init(i,:), & vicen_init=vicen_init(i,:), & trcrn=trcrn(i,1:ntrcr,:), & aice0=aice0(i), & aice =aice(i), & trcr_depend=trcr_depend(1:ntrcr), & trcr_base=trcr_base(1:ntrcr,:), & n_trcr_strata=n_trcr_strata(1:ntrcr), & nt_strata=nt_strata(1:ntrcr,:), & Tf=Tf(i), sss=sss(i), & salinz=salinz(i,:), fside=fside(i), & rside=rside(i), meltl=meltl(i), & frzmlt=frzmlt(i), frazil=frazil(i), & frain=frain(i), fpond=fpond(i), & fresh=fresh(i), fsalt=fsalt(i), & fhocn=fhocn(i), update_ocn_f=update_ocn_f, & bgrid=bgrid, cgrid=cgrid, & igrid=igrid, faero_ocn=faero_ocn(i,:), & first_ice=first_ice(i,:), & fzsal=fzsal(i), & flux_bio=flux_bio(i,1:nbtrcr), & ocean_bio=ocean_bio(i,1:nbtrcr), & frazil_diag=frazil_diag(i), & frz_onset=frz_onset(i), & yday=yday, & nfsd=nfsd, wave_sig_ht=wave_sig_ht(i), & wave_spectrum=wave_spectrum(i,:), & wavefreq=wavefreq(:), & dwavefreq=dwavefreq(:), & d_afsd_latg=d_afsd_latg(i,:), & d_afsd_newi=d_afsd_newi(i,:), & d_afsd_latm=d_afsd_latm(i,:), & d_afsd_weld=d_afsd_weld(i,:), & floe_rad_c=floe_rad_c(:), & floe_binwidth=floe_binwidth(:)) enddo ! i call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) end subroutine step_therm2 !======================================================================= ! ! finalize thermo updates ! subroutine update_state (dt, daidt, dvidt, dagedt, offset) ! column package includes use icepack_intfc, only: icepack_aggregate use icepack_intfc, only: icepack_init_trcr implicit none real (kind=dbl_kind), intent(in) :: & dt , & ! time step offset ! d(age)/dt time offset = dt for thermo, 0 for dyn real (kind=dbl_kind), dimension(:), intent(inout) :: & daidt, & ! change in ice area per time step dvidt, & ! change in ice volume per time step dagedt ! change in ice age per time step integer (kind=int_kind) :: & i, n, k, & ! horizontal indices ntrcr, & ! nt_iage ! real (kind=dbl_kind) :: & temp, & tempsum, & ! tempv, temps, tempa, pi, & Tsfc, sum, hbar, & rhos, Lfresh, puny ! logical (kind=log_kind) :: tr_brine, tr_lvl, tr_fsd integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno, nt_sice, nt_fsd integer (kind=int_kind) :: nt_fbri, nt_alvl, nt_vlvl real (kind=dbl_kind), dimension(ncat) :: & ainit, hinit ! initial area, thickness real (kind=dbl_kind), dimension(nilyr) :: & qin ! ice enthalpy (J/m3) real (kind=dbl_kind), dimension(nslyr) :: & qsn ! snow enthalpy (J/m3) logical (kind=log_kind) :: & tr_iage ! ice age tracer character(len=*), parameter :: subname='(update_state)' !----------------------------------------------------------------- ! query icepack values !----------------------------------------------------------------- call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_lvl_out=tr_lvl, & tr_fsd_out=tr_fsd) call icepack_query_tracer_sizes(ntrcr_out=ntrcr) call icepack_query_tracer_indices( nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, & nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fsd_out=nt_fsd, & nt_fbri_out=nt_fbri, nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl) call icepack_query_parameters(rhos_out=rhos, Lfresh_out=Lfresh, puny_out=puny, pi_out=pi) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_tracer_indices(nt_iage_out=nt_iage) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_tracer_flags(tr_iage_out=tr_iage) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) ! set permanent ice cap at the north pole ainit(:) = c0 hinit(:) = c0 if (3 <= ncat) then n = 3 ainit(n) = c1 ! assumes we are using the default ITD boundaries hinit(n) = c2 else ainit(ncat) = c1 hinit(ncat) = c2 endif ! do i = 1, nx ! do n = 1, ncat ! if((lat_val(i)*180/pi)>91) then ! !write(*,*) i,lat_val(i)*180/pi,lat_val(i) ! aicen(i,n) = ainit(n) ! vicen(i,n) = hinit(n) * ainit(n) ! m ! vsnon(i,n) = c0 ! call icepack_init_trcr(Tair = T_air(i), & ! Tf = Tf(i), & ! Sprofile = salinz(i,:), & ! Tprofile = Tmltz(i,:), & ! Tsfc = Tsfc, & ! nilyr=nilyr, nslyr=nslyr, & ! qin=qin(:), qsn=qsn(:)) ! ! surface temperature ! !trcrn(i,nt_Tsfc,n) = Tsfc ! deg C ! ! ice enthalpy, salinity ! !do k = 1, nilyr ! ! trcrn(i,nt_qice+k-1,n) = qin(k) ! ! trcrn(i,nt_sice+k-1,n) = salinz(i,k) ! !enddo ! ! snow enthalpy ! !do k = 1, nslyr ! ! trcrn(i,nt_qsno+k-1,n) = -rhos * Lfresh ! !enddo ! nslyr ! ! brine fraction ! !if (tr_brine) trcrn(i,nt_fbri,n) = c1 ! endif ! enddo ! ncat ! call icepack_warnings_flush(ice_stderr) ! if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & ! file=__FILE__, line=__LINE__) !enddo !----------------------------------------------------------------- do i = 1, nx !----------------------------------------------------------------- ! Aggregate the updated state variables (includes ghost cells). !----------------------------------------------------------------- call icepack_aggregate (ncat=ncat, & aicen=aicen(i,:), trcrn=trcrn(i,1:ntrcr,:), & vicen=vicen(i,:), vsnon=vsnon(i,:), & aice =aice (i), trcr =trcr (i,1:ntrcr), & vice =vice (i), vsno =vsno (i), & aice0=aice0(i), & ntrcr=ntrcr, & trcr_depend=trcr_depend (1:ntrcr), & trcr_base=trcr_base (1:ntrcr,:), & n_trcr_strata=n_trcr_strata(1:ntrcr), & nt_strata=nt_strata (1:ntrcr,:)) !----------------------------------------------------------------- ! Compute thermodynamic area and volume tendencies. !----------------------------------------------------------------- daidt(i) = (aice(i) - daidt(i)) / dt dvidt(i) = (vice(i) - dvidt(i)) / dt if (tr_iage) then if (offset > c0) then ! thermo if (trcr(i,nt_iage) > c0) & dagedt(i) = (trcr(i,nt_iage) & - dagedt(i) - offset) / dt else ! dynamics dagedt(i) = (trcr(i,nt_iage) & - dagedt(i)) / dt endif endif enddo ! i call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) end subroutine update_state !======================================================================= ! ! Run one time step of wave-fracturing the floe size distribution ! subroutine step_dyn_wave (dt) ! column package includes use icepack_intfc, only: icepack_step_wavefracture use schism_glbl, only : idry implicit none real (kind=dbl_kind), intent(in) :: & dt ! time step ! local variables integer (kind=int_kind) :: & i, j, & ! horizontal indices ntrcr, & ! nbtrcr ! character (len=char_len) :: wave_spec_type character(len=*), parameter :: subname = '(step_dyn_wave)' call icepack_query_parameters(wave_spec_type_out=wave_spec_type) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) do i = 1, nx !if(idry(i)==1) cycle d_afsd_wave(i,:) = c0 call icepack_step_wavefracture (wave_spec_type=wave_spec_type, & dt=dt, ncat=ncat, nfsd=nfsd, nfreq=nfreq, & aice = aice (i), & vice = vice (i), & aicen = aicen (i,:), & floe_rad_l = floe_rad_l (:), & floe_rad_c = floe_rad_c (:), & wave_spectrum = wave_spectrum(i,:), & wavefreq = wavefreq (:), & dwavefreq = dwavefreq (:), & trcrn = trcrn (i,:,:), & d_afsd_wave = d_afsd_wave (i,:)) end do ! i call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) end subroutine step_dyn_wave !======================================================================= ! ! Run one time step of ridging. ! subroutine step_dyn_ridge (dt, ndtd) ! column package includes use icepack_intfc, only: icepack_step_ridge use schism_glbl, only : idry_e,nne,indel implicit none real (kind=dbl_kind), intent(in) :: & dt ! time step integer (kind=int_kind), intent(in) :: & ndtd ! number of dynamics subcycles ! local variables integer (kind=int_kind) :: & i, & ! horizontal indices ntrcr, & ! nbtrcr, & ! narr,n character(len=*), parameter :: subname='(step_dyn_ridge)' !----------------------------------------------------------------- ! query icepack values !----------------------------------------------------------------- call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_out=nbtrcr) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) !----------------------------------------------------------------- ! Ridging !----------------------------------------------------------------- narr = 1 + ncat * (3 + ntrcr) ! max number of state variable arrays do i = 1, nx !if(maxval(idry_e(indel(1:nne(i),i)))/=0) cycle call icepack_step_ridge (dt=dt, ndtd=ndtd, & nilyr=nilyr, nslyr=nslyr, & nblyr=nblyr, & ncat=ncat, hin_max=hin_max(:), & rdg_conv=rdg_conv(i), rdg_shear=rdg_shear(i), & aicen=aicen(i,:), & trcrn=trcrn(i,1:ntrcr,:), & vicen=vicen(i,:), vsnon=vsnon(i,:), & aice0=aice0(i), & trcr_depend=trcr_depend(1:ntrcr), & trcr_base=trcr_base(1:ntrcr,:), & n_trcr_strata=n_trcr_strata(1:ntrcr), & nt_strata=nt_strata(1:ntrcr,:), & dardg1dt=dardg1dt(i), dardg2dt=dardg2dt(i), & dvirdgdt=dvirdgdt(i), opening=opening(i), & fpond=fpond(i), & fresh=fresh(i), fhocn=fhocn(i), & n_aero=n_aero, & faero_ocn=faero_ocn(i,:), & aparticn=aparticn(i,:), krdgn=krdgn(i,:), & aredistn=aredistn(i,:), vredistn=vredistn(i,:), & dardg1ndt=dardg1ndt(i,:), dardg2ndt=dardg2ndt(i,:), & dvirdgndt=dvirdgndt(i,:), & araftn=araftn(i,:), vraftn=vraftn(i,:), & aice=aice(i), fsalt=fsalt(i), & first_ice=first_ice(i,:), fzsal=fzsal(i), & flux_bio=flux_bio(i,1:nbtrcr)) enddo call cut_off_icepack call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) end subroutine step_dyn_ridge !======================================================================= ! ! Computes radiation fields ! subroutine step_radiation (dt) ! column package includes use icepack_intfc, only: icepack_step_radiation use mice_therm_mod, only:t_oi implicit none real (kind=dbl_kind), intent(in) :: & dt ! time step ! local variables integer (kind=int_kind) :: & i, n, k ! horizontal indices integer (kind=int_kind) :: & max_aero, max_algae, nt_Tsfc, nt_alvl, & nt_apnd, nt_hpnd, nt_ipnd, nt_aero, nlt_chl_sw, & ntrcr, nbtrcr_sw, nt_fbri integer (kind=int_kind), dimension(:), allocatable :: & nlt_zaero_sw, nt_zaero, nt_bgc_N logical (kind=log_kind) :: & tr_bgc_N, tr_zaero, tr_brine, dEdd_algae, modal_aero real (kind=dbl_kind), dimension(ncat) :: & fbri ! brine height to ice thickness real(kind= dbl_kind), dimension(:,:), allocatable :: & ztrcr_sw logical (kind=log_kind) :: & l_print_point ! flag for printing debugging information character(len=*), parameter :: subname='(step_radiation)' !----------------------------------------------------------------- ! query icepack values !----------------------------------------------------------------- call icepack_query_tracer_sizes( & max_aero_out=max_aero, max_algae_out=max_algae) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) allocate(nlt_zaero_sw(max_aero)) allocate(nt_zaero(max_aero)) allocate(nt_bgc_N(max_algae)) call icepack_query_tracer_sizes(ntrcr_out=ntrcr, nbtrcr_sw_out=nbtrcr_sw) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_tracer_flags( & tr_brine_out=tr_brine, tr_bgc_N_out=tr_bgc_N, tr_zaero_out=tr_zaero) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_tracer_indices( & nt_Tsfc_out=nt_Tsfc, nt_alvl_out=nt_alvl, nt_apnd_out=nt_apnd, & nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, & nlt_chl_sw_out=nlt_chl_sw, nlt_zaero_sw_out=nlt_zaero_sw, & nt_fbri_out=nt_fbri, nt_zaero_out=nt_zaero, nt_bgc_N_out=nt_bgc_N) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) call icepack_query_parameters(dEdd_algae_out=dEdd_algae, modal_aero_out=modal_aero) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) !----------------------------------------------------------------- allocate(ztrcr_sw(nbtrcr_sw,ncat)) l_print_point = .false. do i = 1, nx fbri(:) = c0 ztrcr_sw(:,:) = c0 do n = 1, ncat if (tr_brine) fbri(n) = trcrn(i,nt_fbri,n) enddo call icepack_step_radiation(dt=dt, ncat=ncat, & nblyr=nblyr, nilyr=nilyr, & nslyr=nslyr, dEdd_algae=dEdd_algae, & swgrid=swgrid(:), igrid=igrid(:), & fbri=fbri(:), & aicen=aicen(i,:), vicen=vicen(i,:), & vsnon=vsnon(i,:), & Tsfcn=trcrn(i,nt_Tsfc,:), & alvln=trcrn(i,nt_alvl,:), & apndn=trcrn(i,nt_apnd,:), & hpndn=trcrn(i,nt_hpnd,:), & ipndn=trcrn(i,nt_ipnd,:), & aeron=trcrn(i,nt_aero:nt_aero+4*n_aero-1,:), & bgcNn=trcrn(i,nt_bgc_N(1):nt_bgc_N(1)+n_algae*(nblyr+3)-1,:), & zaeron=trcrn(i,nt_zaero(1):nt_zaero(1)+n_zaero*(nblyr+3)-1,:), & trcrn_bgcsw=ztrcr_sw, & TLAT=lat_val(i), TLON=lon_val(i), & calendar_type=calendar_type, & days_per_year=days_per_year, sec=sec, & nextsw_cday=nextsw_cday, yday=yday, & kaer_tab=kaer_tab, kaer_bc_tab=kaer_bc_tab(:,:), & waer_tab=waer_tab, waer_bc_tab=waer_bc_tab(:,:), & gaer_tab=gaer_tab, gaer_bc_tab=gaer_bc_tab(:,:), & bcenh=bcenh(:,:,:), modal_aero=modal_aero, & swvdr=swvdr(i), swvdf=swvdf(i), & swidr=swidr(i), swidf=swidf(i), & coszen=cos_zen(i), fsnow=fsnow(i), & alvdrn=alvdrn(i,:), alvdfn=alvdfn(i,:), & alidrn=alidrn(i,:), alidfn=alidfn(i,:), & fswsfcn=fswsfcn(i,:), fswintn=fswintn(i,:), & fswthrun=fswthrun(i,:), fswpenln=fswpenln(i,:,:), & Sswabsn=Sswabsn(i,:,:), Iswabsn=Iswabsn(i,:,:), & albicen=albicen(i,:), albsnon=albsnon(i,:), & albpndn=albpndn(i,:), apeffn=apeffn(i,:), & snowfracn=snowfracn(i,:), & dhsn=dhsn(i,:), ffracn=ffracn(i,:), & l_print_point=l_print_point) if (dEdd_algae .and. (tr_zaero .or. tr_bgc_N)) then do n = 1, ncat do k = 1, nbtrcr_sw trcrn_sw(i,k,n) = ztrcr_sw(k,n) enddo enddo endif t_oi(i)=trcr(i,nt_Tsfc)!t_oi enddo ! i call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) deallocate(ztrcr_sw) deallocate(nlt_zaero_sw) deallocate(nt_zaero) deallocate(nt_bgc_N) end subroutine step_radiation !======================================================================= ! ! Ocean mixed layer calculation (internal to sea ice model). ! Allows heat storage in ocean for uncoupled runs. ! subroutine ocean_mixed_layer (dt) use icepack_intfc, only: icepack_ocn_mixed_layer, icepack_atm_boundary implicit none real (kind=dbl_kind), intent(in) :: & dt ! time step ! local variables integer (kind=int_kind) :: & i ! horizontal indices real (kind=dbl_kind) :: & albocn real (kind=dbl_kind), dimension(nx) :: & delt , & ! potential temperature difference (K) delq , & ! specific humidity difference (kg/kg) shcoef, & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat character(len=*), parameter :: subname='(ocean_mixed_layer)' !----------------------------------------------------------------- ! query icepack values !----------------------------------------------------------------- call icepack_query_parameters(albocn_out=albocn) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) !----------------------------------------------------------------- ! Compute boundary layer quantities !----------------------------------------------------------------- do i = 1, nx call icepack_atm_boundary(sfctype = 'ocn', & Tsf = sst(i), & potT = potT(i), & uatm = uatm(i), & vatm = vatm(i), & wind = wind(i), & !zlvl_t = zlvl_t, & !zlvl_q = zlvl_q, & zlvl = zlvl_v, & Qa = Qa(i), & rhoa = rhoa(i), & strx = strairx_ocn(i), & stry = strairy_ocn(i), & Tref = Tref_ocn(i), & Qref = Qref_ocn(i), & delt = delt(i), & delq = delq(i), & lhcoef = lhcoef(i), & shcoef = shcoef(i), & Cdn_atm = Cdn_atm(i), & Cdn_atm_ratio_n = Cdn_atm_ratio(i)) enddo ! i call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) !----------------------------------------------------------------- ! Ocean albedo ! For now, assume albedo = albocn in each spectral band. !----------------------------------------------------------------- alvdr_ocn(:) = albocn alidr_ocn(:) = albocn alvdf_ocn(:) = albocn alidf_ocn(:) = albocn !----------------------------------------------------------------- ! Compute ocean fluxes and update SST !----------------------------------------------------------------- do i = 1, nx call ocn_mixed_layer_icepack( alvdr_ocn=alvdr_ocn(i), swvdr=swvdr(i), & alidr_ocn=alidr_ocn(i), swidr=swidr(i), & alvdf_ocn=alvdf_ocn(i), swvdf=swvdf(i), & alidf_ocn=alidf_ocn(i), swidf=swidf(i), & flwout_ocn=flwout_ocn(i),sst=sst(i), & fsens_ocn=fsens_ocn(i), shcoef=shcoef(i), & flat_ocn=flat_ocn(i), lhcoef=lhcoef(i), & evap_ocn=evap_ocn(i), flw=flw(i), & delt=delt(i), delq=delq(i), & aice=aice(i), fhocn=fhocn(i), & fswthru=fswthru(i), hmix=hmix(i), & Tf=Tf(i), fresh=fresh(i), & frain=frain(i), fsnow=fsnow(i), & fhocn_tot=fhocn_tot(i), fresh_tot=fresh_tot(i), & frzmlt=frzmlt(i), fsalt=fsalt(i), & sss=sss(i), dt =dt ) enddo ! i call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) end subroutine ocean_mixed_layer !======================================================================= subroutine ocn_mixed_layer_icepack( & alvdr_ocn, swvdr, & alidr_ocn, swidr, & alvdf_ocn, swvdf, & alidf_ocn, swidf, & sst, flwout_ocn, & fsens_ocn, shcoef, & flat_ocn, lhcoef, & evap_ocn, flw, & delt, delq, & aice, fhocn, & fswthru, hmix, & Tf, fresh, & frain, fsnow, & fhocn_tot, fresh_tot, & frzmlt, fsalt, & sss, dt) !use i_therm_param, only: emiss_wat use mice_therm_mod, only: emiss_wat !use g_forcing_param, only: use_virt_salt implicit none logical :: use_virt_salt=.false.! will be set TRUE in case of which_ALE='linfs', otherwise FALSE real (kind=dbl_kind), intent(in) :: & alvdr_ocn , & ! visible, direct (fraction) alidr_ocn , & ! near-ir, direct (fraction) alvdf_ocn , & ! visible, diffuse (fraction) alidf_ocn , & ! near-ir, diffuse (fraction) swvdr , & ! sw down, visible, direct (W/m^2) swvdf , & ! sw down, visible, diffuse (W/m^2) swidr , & ! sw down, near IR, direct (W/m^2) swidf , & ! sw down, near IR, diffuse (W/m^2) flw , & ! incoming longwave radiation (W/m^2) Tf , & ! freezing temperature (C) hmix , & ! mixed layer depth (m) delt , & ! potential temperature difference (K) delq , & ! specific humidity difference (kg/kg) shcoef , & ! transfer coefficient for sensible heat lhcoef , & ! transfer coefficient for latent heat fswthru , & ! shortwave penetrating to ocean (W/m^2) aice , & ! ice area fraction sst , & ! sea surface temperature (C) sss , & ! sea surface salinity frain , & ! rainfall rate (kg/m^2/s) fsnow , & ! snowfall rate (kg/m^2/s) fsalt , & ! salt flux from ice to the ocean (kg/m^2/s) dt real (kind=dbl_kind), intent(inout) :: & flwout_ocn, & ! outgoing longwave radiation (W/m^2) fsens_ocn , & ! sensible heat flux (W/m^2) flat_ocn , & ! latent heat flux (W/m^2) evap_ocn , & ! evaporative water flux (kg/m^2/s) fhocn , & ! net heat flux to ocean (W/m^2) fresh , & ! fresh water flux to ocean (kg/m^2/s) frzmlt , & ! freezing/melting potential (W/m^2) fhocn_tot , & ! net total heat flux to ocean (W/m^2) fresh_tot ! fresh total water flux to ocean (kg/m^2/s) real (kind=dbl_kind), parameter :: & frzmlt_max = c1000 ! max magnitude of frzmlt (W/m^2) real (kind=dbl_kind) :: & TsfK , & ! surface temperature (K) swabs, & ! surface absorbed shortwave heat flux (W/m^2) Tffresh, & ! 0 C in K Lfresh, & Lvap, & lfs_corr, & ! fresh water correction for linear free surface stefan_boltzmann, & ice_ref_salinity, & cprho character(len=*),parameter :: subname='(icepack_ocn_mixed_layer)' call icepack_query_parameters( Tffresh_out=Tffresh, Lfresh_out=Lfresh, & stefan_boltzmann_out=stefan_boltzmann, & ice_ref_salinity_out=ice_ref_salinity, & Lvap_out=Lvap,cprho_out=cprho ) ! shortwave radiative flux ! Visible is absorbed by clorophil ! afterwards swabs = (c1-alidr_ocn) * swidr + (c1-alidf_ocn) * swidf + & (c1-alvdr_ocn) * swvdr + (c1-alvdf_ocn) * swvdf ! ocean surface temperature in Kelvin TsfK = sst + Tffresh ! longwave radiative flux ! Water emissivity added to be consistent ! with the standard FESOM2 version flwout_ocn = - emiss_wat * stefan_boltzmann * TsfK**4 ! downward latent and sensible heat fluxes fsens_ocn = shcoef * delt flat_ocn = lhcoef * delq evap_ocn = -flat_ocn / Lvap ! Compute heat change due to exchange between ocean and atmosphere fhocn_tot = fhocn & ! these are *aice already + (fsens_ocn + flat_ocn + flwout_ocn + flw & - Lfresh*fsnow) * (c1-aice)+ max(c0,frzmlt)!*aice !- Lfresh*fsnow) * (c1-aice)+ max(c0,frzmlt)!*aice !+ swabs & !+ Lfresh*fsnow) * (c1-aice) !+ max(c0,frzmlt)*aice !+ min(max(-frzmlt_max,frzmlt),frzmlt_max)*aice! !fhocn_tot = fhocn + fswthru & ! these are *aice already ! + (fsens_ocn + flat_ocn + flwout_ocn + flw + swabs & ! + Lfresh*fsnow) * (c1-aice) + max(c0,frzmlt)*aice if (use_virt_salt) then lfs_corr = fsalt/ice_ref_salinity/p001 fresh = fresh - lfs_corr * ice_ref_salinity / sss endif fresh = fresh - max(fsalt,0.d0) /p001 / sss fresh_tot = fresh + (-evap_ocn + frain + fsnow)*(c1-aice) end subroutine ocn_mixed_layer_icepack !======================================================================= subroutine coupling_prep(dt) ! local variables implicit none real (kind=dbl_kind), intent(in) :: & dt ! time step integer (kind=int_kind) :: & n , & ! thickness category index i , & ! horizontal index k , & ! tracer index nbtrcr real (kind=dbl_kind) :: & netsw, & ! flag for shortwave radiation presence rhofresh, & ! puny ! character(len=*), parameter :: subname='(coupling_prep)' !----------------------------------------------------------------- ! Save current value of frzmlt for diagnostics. ! Update mixed layer with heat and radiation from ice. !----------------------------------------------------------------- call icepack_query_parameters(puny_out=puny, rhofresh_out=rhofresh) call icepack_query_tracer_sizes(nbtrcr_out=nbtrcr) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) do i = 1, nx frzmlt_init (i) = frzmlt(i) enddo call ocean_mixed_layer (dt) ! ocean surface fluxes and sst !----------------------------------------------------------------- ! Aggregate albedos !----------------------------------------------------------------- do i = 1, nx alvdf(i) = c0 alidf(i) = c0 alvdr(i) = c0 alidr(i) = c0 albice(i) = c0 albsno(i) = c0 albpnd(i) = c0 apeff_ai(i) = c0 snowfrac(i) = c0 enddo do n = 1, ncat do i = 1, nx if (aicen(i,n) > puny) then alvdf(i) = alvdf(i) + alvdfn(i,n)*aicen(i,n) alidf(i) = alidf(i) + alidfn(i,n)*aicen(i,n) alvdr(i) = alvdr(i) + alvdrn(i,n)*aicen(i,n) alidr(i) = alidr(i) + alidrn(i,n)*aicen(i,n) netsw = swvdr(i) + swidr(i) + swvdf(i) + swidf(i) if (netsw > puny) then ! sun above horizon albice(i) = albice(i) + albicen(i,n)*aicen(i,n) albsno(i) = albsno(i) + albsnon(i,n)*aicen(i,n) albpnd(i) = albpnd(i) + albpndn(i,n)*aicen(i,n) endif apeff_ai(i) = apeff_ai(i) + apeffn(i,n)*aicen(i,n) ! for history snowfrac(i) = snowfrac(i) + snowfracn(i,n)*aicen(i,n) ! for history endif ! aicen > puny enddo enddo do i = 1, nx !----------------------------------------------------------------- ! reduce fresh by fpond for coupling !----------------------------------------------------------------- if (l_mpond_fresh) then fpond(i) = fpond(i) * rhofresh/dt fresh(i) = fresh(i) - fpond(i) endif !---------------------------------------------------------------- ! Store grid box mean albedos and fluxes before scaling by aice !---------------------------------------------------------------- alvdf_ai (i) = alvdf (i) alidf_ai (i) = alidf (i) alvdr_ai (i) = alvdr (i) alidr_ai (i) = alidr (i) fresh_ai (i) = fresh (i) fsalt_ai (i) = fsalt (i) fhocn_ai (i) = fhocn (i) fswthru_ai(i) = fswthru(i) fzsal_ai (i) = fzsal (i) fzsal_g_ai(i) = fzsal_g(i) if (nbtrcr > 0) then do k = 1, nbtrcr flux_bio_ai (i,k) = flux_bio (i,k) enddo endif !----------------------------------------------------------------- ! Save net shortwave for scaling factor in scale_factor !----------------------------------------------------------------- scale_factor(i) = & swvdr(i)*(c1 - alvdr_ai(i)) & + swvdf(i)*(c1 - alvdf_ai(i)) & + swidr(i)*(c1 - alidr_ai(i)) & + swidf(i)*(c1 - alidf_ai(i)) enddo end subroutine coupling_prep !======================================================================= module subroutine step_icepack() use schism_glbl, only: rkind,pi,np,npa,nvrt,uu2,vv2,time_stamp,windx,windy,xnd,ynd, & &tau_oi,nws,ihconsv,isconsv,iplg,fresh_wa_flux,net_heat_flux,srad_th_ice,rho0,rnday,fdb,lfdb, & &lice_free_gb,lhas_ice,errmsg,ice_evap,isbnd,nnp,indnd,nstep_ice,it_main use schism_msgp, only : myrank,nproc,parallel_abort,comm,ierr,exchange_p2d use mice_module use mice_therm_mod implicit none integer (kind=int_kind) :: & i,k,jj,njj,kk, nt_Tsfc,nt_sice ! dynamics supercycling index logical (kind=log_kind) :: & calc_Tsfc, skl_bgc, solve_zsal, z_tracers, tr_brine, & ! from icepack tr_fsd, wave_spec real (kind=dbl_kind) :: & offset, & ! d(age)/dt time offset t1, t2, t3, t4,dt,umod,tmp1,& dux,dvy,aux,puny real(kind=dbl_kind) :: & swild(npa) logical :: lice_free !real (kind=dbl_kind), intent(out) :: & ! time_therm, & ! time_advec, & ! time_evp !type(t_mesh), target, intent(in) :: mesh character(len=*), parameter :: subname='(ice_step)' dt=ice_dt !=dt*nstep_ice !t1 = c0 !t2 = c0 !t3 = c0 !t4 = c0 !t1 = MPI_Wtime() !----------------------------------------------------------------- ! query Icepack values !----------------------------------------------------------------- call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers) call icepack_query_parameters(solve_zsal_out=solve_zsal, calc_Tsfc_out=calc_Tsfc, & wave_spec_out=wave_spec,puny_out=puny) call icepack_query_tracer_indices( & nt_Tsfc_out=nt_Tsfc,nt_sice_out=nt_sice) call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_fsd_out=tr_fsd) call icepack_warnings_flush(ice_stderr) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) ! TODO: Add appropriate timing !----------------------------------------------------------------- ! copy variables from schism (also ice velocities) !----------------------------------------------------------------- call schism_to_icepack !----------------------------------------------------------------- ! tendencies needed by schism !----------------------------------------------------------------- dhi_dt(:) = vice(:) dhs_dt(:) = vsno(:) !----------------------------------------------------------------- ! initialize diagnostics !----------------------------------------------------------------- call init_history_therm call init_history_bgc !----------------------------------------------------------------- ! Scale radiation fields !----------------------------------------------------------------- if (calc_Tsfc) call prep_radiation () !----------------------------------------------------------------- ! thermodynamics and biogeochemistry !----------------------------------------------------------------- if(mod(it_main-1,nstep_ice)==0) then if(ice_therm_on==1) then if(ice_tests==0.and.(nws/=2.or.ihconsv/=1.or.isconsv/=1)) & &call parallel_abort('ice_step: ice therm needs nws=2 etc') !Atmos variables are read in for thermodynamics call step_therm1 (dt) ! vertical thermodynamics call step_therm2 (dt) ! ice thickness distribution thermo if(myrank==0) write(16,*)'done ice thermodynamics' endif endif ! clean up, update tendency diagnostics offset = dt call update_state (dt, daidtt, dvidtt, dagedtt, offset) !----------------------------------------------------------------- ! dynamics, transport, ridging !----------------------------------------------------------------- call init_history_dyn ! wave fracture of the floe size distribution ! note this is called outside of the dynamics subcycling loop if (tr_fsd .and. wave_spec) call step_dyn_wave(dt) do k = 1, ndtd !split in time; =1 in icedrv_set.F90 !----------------------------------------------------------------- ! EVP !----------------------------------------------------------------- !t2 = MPI_Wtime() if(ievp==1) then call ice_evp if(myrank==0) write(16,*)'done ice EVP dynamics' else if (ievp==2) then call ice_mevp if(myrank==0) write(16,*)'done ice mEVP dynamics' else if(myrank==0) write(16,*)'no ice dynamics' endif !ievp !Transport: operator splitting !if(ice_advection/=0) then ! call ice_fct ! if(myrank==0) write(16,*)'done ice FCT advection' !endif !t3 = MPI_Wtime() !time_evp = t3 - t2 lice_free=.true. !over all aug domain call icepack_to_schism (nx_in=npa, & aice_out=a_ice0, & vice_out=m_ice0, & vsno_out=m_snow0, & fhocn_tot_out=net_heat_flux0, & fresh_tot_out=fresh_wa_flux0, & strocnxT_out=tau_oi_x1, & strocnyT_out=tau_oi_y1, & !fsalt_out=real_salt_flux, & !dhi_dt_out=thdgrsn, & !dhs_dt_out=thdgr, & evap_ocn_out=evaporation ) do i=1,npa if(a_ice0(i)<=ice_cutoff.or.m_ice0(i)<=ice_cutoff) then lhas_ice(i)=.false. U_ice(i)=0 V_ice(i)=0 else lhas_ice(i)=.true. lice_free=.false. endif enddo call mpi_allreduce(lice_free,lice_free_gb,1,MPI_LOGICAL,MPI_LAND,comm,ierr) if(myrank==0) write(16,*)'lice_free_gb=',lice_free_gb !----------------------------------------------------------------- ! update ice velocities !----------------------------------------------------------------- call schism_to_icepack !m_ice=ice_tr(1,:) !a_ice=ice_tr(2,:) !m_snow=ice_tr(3,:) !----------------------------------------------------------------- ! advect tracers !----------------------------------------------------------------- !t2 = MPI_Wtime() if(ice_advection==1) then call tracer_advection_icepack if(myrank==0) write(16,*)'done ice FCT advection' else if (ice_advection==2) then call tracer_advection_icepack2 if(myrank==0) write(16,*)'done ice advection in Taylor-Galerkin way' else if (ice_advection==3) then call tracer_advection_icepack3 if(myrank==0) write(16,*)'done ice upwind advection' else if(myrank==0) write(16,*)'no ice dynamics' endif !t3 = MPI_Wtime() !time_advec = t3 - t2 !----------------------------------------------------------------- ! ridging !----------------------------------------------------------------- call step_dyn_ridge (dt_dyn, ndtd) ! clean up, update tendency diagnostics offset = c0 call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset) enddo !call schism_to_icepack !----------------------------------------------------------------- ! albedo, shortwave radiation !----------------------------------------------------------------- call step_radiation (dt) !----------------------------------------------------------------- ! get ready for coupling and the next time step !----------------------------------------------------------------- call coupling_prep (dt) !----------------------------------------------------------------- ! tendencies needed by fesom !----------------------------------------------------------------- dhi_dt(:) = ( vice(:) - dhi_dt(:) ) / dt dhs_dt(:) = ( vsno(:) - dhi_dt(:) ) / dt !----------------------------------------------------------------- ! icepack timing !----------------------------------------------------------------- !t4 = MPI_Wtime() !time_therm = t4 - t1 - time_advec - time_evp !time_advec = c0 !time_therm = c0 !time_evp = c0 call icepack_to_schism (nx_in=npa, & aice_out=a_ice0, & vice_out=m_ice0, & vsno_out=m_snow0, & fhocn_tot_out=net_heat_flux0, & fresh_tot_out=fresh_wa_flux0, & strocnxT_out=tau_oi_x1, & strocnyT_out=tau_oi_y1, & !fsalt_out=real_salt_flux, & !dhi_dt_out=thdgrsn, & !dhs_dt_out=thdgr, & evap_ocn_out=evaporation, & fsrad_ice_out=fsrad_ice_out0 ) !Mark ice/no ice lice_free=.true. !over all aug domain do i=1,npa if(a_ice0(i)<=ice_cutoff.or.m_ice0(i)<=ice_cutoff) then lhas_ice(i)=.false. U_ice(i)=0 V_ice(i)=0 else lhas_ice(i)=.true. lice_free=.false. endif !write(12,*) 'ice',i,ice_cutoff,a_ice0(i),m_ice0(i),m_snow0(i),lhas_ice(i) enddo !i call mpi_allreduce(lice_free,lice_free_gb,1,MPI_LOGICAL,MPI_LAND,comm,ierr) if(myrank==0) write(16,*)'lice_free_gb=',lice_free_gb fresh_wa_flux(:)=0 net_heat_flux(:)=0 tau_oi=0 !init as junk ice_tr(:,:)=0 ice_evap(:)=0 srad_th_ice(:)=0 do i=1,npa !if(lhas_ice(i)) then !tau_oi(1,i)=tau_oi_x1(i) !m^2/s/s !tau_oi(2,i)=tau_oi_y1(i) ice_tr(1,i)=m_ice0(i) ice_tr(2,i)=a_ice0(i) ice_tr(3,i)=m_snow0(i) net_heat_flux(i) = net_heat_flux0(i) fresh_wa_flux(i) = fresh_wa_flux0(i) !- runoff(:) ice_evap(i)=evaporation(i)*(1-a_ice0(i)) srad_th_ice(i) = fsrad_ice_out0(i) !endif !Check outputs !if(t_oi(i)/=t_oi(i).or.fresh_wa_flux(i)/=fresh_wa_flux(i).or. & !&net_heat_flux(i)/=net_heat_flux(i).or.ice_tr(1,i)/=ice_tr(1,i).or. & !&ice_tr(2,i)/=ice_tr(2,i).or.ice_tr(3,i)/=ice_tr(3,i).or.ice_evap(i)/=ice_evap(i)) then ! write(errmsg,*)'ice_thermodynamics: nan',iplg(i),t_oi(i), & !&fresh_wa_flux(i),net_heat_flux(i),ice_tr(:,i) ! call parallel_abort(errmsg) ! endif enddo !i=1,npa do i=1,np if(lhas_ice(i)) then umod=sqrt((U_ice(i)-u_ocean(i))**2+(V_ice(i)-v_ocean(i))**2) tmp1=ice_tr(2,i)*Cdn_ocn(i)*umod if(isbnd(1,i)<0) then !b.c. (including open) njj=0 do kk=1,nnp(i) jj=indnd(kk,i) if((isbnd(1,jj)==0).and.(lhas_ice(jj))) then umod=sqrt((U_ice(jj)-u_ocean(jj))**2+(V_ice(jj)-v_ocean(jj))**2) tmp1=a_ice0(jj)*Cdn_ocn(i)*umod tau_oi(1,i)=tau_oi(1,i)+tmp1*((U_ice(jj)-u_ocean(jj))*cos_io-(V_ice(jj)-v_ocean(jj))*sin_io) !m^2/s/s tau_oi(2,i)=tau_oi(2,i)+tmp1*((U_ice(jj)-u_ocean(jj))*sin_io+(V_ice(jj)-v_ocean(jj))*cos_io) njj=njj+1 endif enddo tau_oi(1,i)=tau_oi(1,i)/njj tau_oi(2,i)=tau_oi(2,i)/njj !tau_oi(1,i)=0.01*tmp1*((U_ice(i)-u_ocean(i))*cos_io-(V_ice(i)-v_ocean(i))*sin_io) !m^2/s/s !tau_oi(2,i)=0.01*tmp1*((U_ice(i)-u_ocean(i))*sin_io+(V_ice(i)-v_ocean(i))*cos_io) !tau_oi(1,i)=0 !0.001*((U_ice(i)-u_ocean(i))*cos_io-(V_ice(i)-v_ocean(i))*sin_io)/umod !tau_oi(2,i)=0 !0.001*((U_ice(i)-u_ocean(i))*sin_io+(V_ice(i)-v_ocean(i))*cos_io)/umod else if(isbnd(1,i)>0) then tau_oi(1,i)=0 tau_oi(2,i)=0 else tau_oi(1,i)=tmp1*((U_ice(i)-u_ocean(i))*cos_io-(V_ice(i)-v_ocean(i))*sin_io) !m^2/s/s tau_oi(2,i)=tmp1*((U_ice(i)-u_ocean(i))*sin_io+(V_ice(i)-v_ocean(i))*cos_io) endif endif if((tau_oi(1,i)/=tau_oi(1,i)).or.(tau_oi(2,i)/=tau_oi(2,i))) then tau_oi(1,i)=0 tau_oi(2,i)=0 endif enddo swild=tau_oi(1,:) call exchange_p2d(swild) tau_oi(1,:)=swild swild=tau_oi(2,:) call exchange_p2d(swild) tau_oi(2,:)=swild !call exchange_p2d(net_heat_flux) !call exchange_p2d(fresh_wa_flux) !call exchange_p2d(ice_evap) !swild=ice_tr(1,:) !call exchange_p2d(swild) !ice_tr(1,:)=swild !swild=ice_tr(2,:) !call exchange_p2d(swild) !ice_tr(2,:)=swild !swild=ice_tr(3,:) !call exchange_p2d(swild) !ice_tr(3,:)=swild call init_flux_atm_ocn() end subroutine step_icepack !======================================================================= end submodule icedrv_step !=======================================================================