!/===========================================================================/ ! Copyright (c) 2007, The University of Massachusetts Dartmouth ! Produced at the School of Marine Science & Technology ! Marine Ecosystem Dynamics Modeling group ! All rights reserved. ! ! FVCOM has been developed by the joint UMASSD-WHOI research team. For ! details of authorship and attribution of credit please see the FVCOM ! technical manual or contact the MEDM group. ! ! ! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu ! The full copyright notice is contained in the file COPYRIGHT located in the ! root directory of the FVCOM code. This original header must be maintained ! in all distributed versions. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, ! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR ! PURPOSE ARE DISCLAIMED. ! !/---------------------------------------------------------------------------/ ! CVS VERSION INFORMATION ! $Id$ ! $Name$ ! $Revision$ !/===========================================================================/ ! VERSIONS. Finite volume ! ! Comments: Sea Ice Module ! connect the thermodynamic module(CICE) and Finite-volume dynamics module(EVP) ! ! Current FVCOM dependency ! archive.F: - calls archive_ice ! mod_ncdio.F: - netcdf output includes some variables ! us_fvcom.F: - main calls ice setup !======================================================================= MODULE MOD_ICE # if defined (ICE) USE ALL_VARS USE MOD_PAR USE MOD_UTILS use mod_ice2d ! afm 20180723 -- begin # if defined (HEATING_CALCULATED) USE MOD_HEATFLUX, ONLY : HEATING_FRESHWATER # endif ! --end !==================================================================== !-------------------------------------------------------------- ! for cice code ! !USES: ! use ice_model_size use ice_albedo use ice_calendar ! use ice_coupling ! use ice_diagnostics use ice_domain ! use ice_dyn_evp use ice_fileunits use ice_flux_in use ice_grid ! use ice_history use ice_init use ice_itd use ice_itd_linear use ice_kinds_mod use ice_mechred # if !defined (ICE_FRESHWATER) ! afm 20160513 & EJA 20160921 not used for freshwater ! use ice_ocean # endif use ice_atmo use ice_scaling use ice_therm_vertical use ice_therm_itd use ice_work !==================================================================== IMPLICIT NONE CONTAINS !===============================================================================| SUBROUTINE ice_init_0 ! !============================================================== IMPLICIT NONE integer :: I,J,K !-------------------------------------------------------------- ! for cice code ! !USES: !======================================================================= !=======================================================================| ! ice_domain_def IF(DBG_SET(DBG_SBR)) write(ipt,*)'START ice_init_0' ilo = 1 ! beg index of actual physical subdomain ihi = M !T ! end index imt_local = M !T !, & jlo = 1 ! beg index of actual physical subdomain jhi =1 jmt_local= 1 allocate(index_global(ilo:ihi,jlo:jhi), rndex_global(ilo:ihi,jlo:jhi)) ! ice_work allocate( work_l1(imt_local,jmt_local)) allocate( work_l2(imt_local,jmt_local)) allocate( worka(ilo:ihi,jlo:jhi)) allocate( workb(ilo:ihi,jlo:jhi)) ! ice_state_def allocate(aice(imt_local,jmt_local)) ;aice=ZERO ! concentration of ice allocate(vice(imt_local,jmt_local)) ;vice=ZERO ! volume per unit area of ice (m) allocate(vsno(imt_local,jmt_local)) ;vsno=ZERO ! volume per unit area of snow (m) allocate(Tsfc(imt_local,jmt_local)) ;Tsfc=ZERO ! temperature of ice/snow top surface (C) allocate(eice(imt_local,jmt_local)) ;eice=ZERO ! energy of melt. of ice (J/m^2) allocate(esno(imt_local,jmt_local)) ;esno=ZERO ! energy of melt. of snow layer (J/m^2) allocate(tmass(imt_local,jmt_local)) ;tmass=ZERO ! total mass of ice and snow allocate(aice_init(imt_local,jmt_local));aice_init=ZERO ! concentration of ice at beginning of dt (for diagnostics) ! - begin afm 20171113 for ice open boundary ------------- allocate(xflux_ice(0:MT));xflux_ice=ZERO ! - end afm 20171113 for ice open boundary ------------- !----------------------------------------------------------------- ! state of the ice for each category !----------------------------------------------------------------- allocate(aicen(imt_local,jmt_local,ncat)) ;aicen=ZERO ! concentration of ice allocate(vicen(imt_local,jmt_local,ncat)) ;vicen=ZERO ! volume per unit area of ice (m) allocate(vsnon(imt_local,jmt_local,ncat)) ;vsnon=ZERO ! volume per unit area of snow (m) allocate(Tsfcn(imt_local,jmt_local,ncat)) ;Tsfcn=ZERO ! temperature of ice/snow top surface (C) allocate(esnon(imt_local,jmt_local,ncat)) ;esnon=ZERO ! energy of melt. of snow layer (J/m^2) allocate(aice0(imt_local,jmt_local)) ;aice0=C1i ! concentration of open water allocate(eicen(imt_local,jmt_local,ntilay));eicen=ZERO ! energy of melting for !----------------------------------------------------------------- ! other variables closely related to the state of the ice !----------------------------------------------------------------- !allocate(uvel(imt_local,jmt_local)) ;uvel=ZERO! x-component of velocity (m/s) !allocate(vvel(imt_local,jmt_local)) ;vvel=ZERO! y-component of velocity (m/s) allocate(strength(imt_local,jmt_local)) ;strength=ZERO ! ice strength (N/m) !======================================================================= ! subroutine ice_albedo_def allocate(alvdrn (ilo:ihi,jlo:jhi,ncat)) ;alvdrn=ZERO ! visible, direct (fraction) allocate(alidrn (ilo:ihi,jlo:jhi,ncat)) ;alidrn=ZERO ! near-ir, direct (fraction) allocate(alvdfn (ilo:ihi,jlo:jhi,ncat)) ;alvdfn=ZERO ! visible, diffuse (fraction) allocate(alidfn (ilo:ihi,jlo:jhi,ncat)) ;alidfn=ZERO ! near-ir, diffuse (fraction) ! albedos aggregated over categories allocate(alvdr (ilo:ihi,jlo:jhi)) ;alvdr=ZERO ! visible, direct (fraction) allocate(alidr (ilo:ihi,jlo:jhi)) ;alidr=ZERO ! near-ir, direct (fraction) allocate(alvdf (ilo:ihi,jlo:jhi)) ;alvdf=ZERO ! visible, diffuse (fraction) allocate(alidf (ilo:ihi,jlo:jhi)) ;alidf=ZERO ! near-ir, diffuse (fraction) !======================================================================= ! subroutine ice_flux_def !----------------------------------------------------------------- ! Dynamics component !----------------------------------------------------------------- ! Dynamic component use FVCOM code ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: & ! in from ocean ! allocate(ss_tltx(ilo:ihi,jlo:jhi)) ;ss_tltx=ZERO ! sea surface slope, x-direction (m/m) ! allocate(ss_tlty(ilo:ihi,jlo:jhi)) ;ss_tlty=ZERO ! sea surface slope, y-direction ! allocate(uocn (ilo:ihi,jlo:jhi)) ;uocn =ZERO ! ocean current, x-direction (m/s) ! allocate(vocn (ilo:ihi,jlo:jhi)) ;vocn =ZERO ! ocean current, y-direction (m/s) ! out to atmosphere allocate(strairxT(ilo:ihi,jlo:jhi)) ;strairxT=ZERO ! stress on ice by air, x-direction allocate(strairyT(ilo:ihi,jlo:jhi)) ;strairyT=ZERO ! stress on ice by air, y-direction ! out to ocean T-cell (kg/m s^2) allocate(strocnxT(ilo:ihi,jlo:jhi)) ;strocnxT=ZERO ! ice-ocean stress, x-direction allocate(strocnyT(ilo:ihi,jlo:jhi)) ;strocnyT=ZERO ! ice-ocean stress, y-direction ! diagnostic allocate(daidtd(ilo:ihi,jlo:jhi)) ;daidtd=ZERO ! ice area tendency due to transport (s^-1) allocate(dvidtd(ilo:ihi,jlo:jhi)) ;dvidtd=ZERO ! ice volume tendency due to transport (m/s) ! internal U-cell (kg/m s^2) allocate(strocnx(ilo:ihi,jlo:jhi)) ;strocnx=ZERO ! ice-ocean stress, x-direction allocate(strocny(ilo:ihi,jlo:jhi)) ;strocny=ZERO ! ice-ocean stress, y-direction allocate(strairx(ilo:ihi,jlo:jhi)) ;strairx=ZERO ! stress on ice by air, x-direction allocate(strairy(ilo:ihi,jlo:jhi)) ;strairy=ZERO ! stress on ice by air, y-direction ! allocate(strtltx(ilo:ihi,jlo:jhi)) ;strtltx=ZERO ! stress due to sea surface slope, x-direction ! allocate(strtlty(ilo:ihi,jlo:jhi)) ;strtlty=ZERO ! stress due to sea surface slope, y-direction !----------------------------------------------------------------- ! Thermodynamics component !----------------------------------------------------------------- ! in from atmosphere ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: & allocate(zlvl (ilo:ihi,jlo:jhi)) ;zlvl =ZERO ! atm level height (m) allocate(uatm (ilo:ihi,jlo:jhi)) ;uatm =ZERO ! wind speed (m/s) allocate(vatm (ilo:ihi,jlo:jhi)) ;vatm =ZERO ! allocate(potT (ilo:ihi,jlo:jhi)) ;potT =ZERO ! air potential temperature (K) allocate(Tair (ilo:ihi,jlo:jhi)) ;Tair =ZERO ! air temperature (K) allocate(Qa (ilo:ihi,jlo:jhi)) ;Qa =ZERO ! specific humidity (kg/kg) allocate(rhoa (ilo:ihi,jlo:jhi)) ;rhoa =1.3_SP! air density (kg/m^3) allocate(swvdr(ilo:ihi,jlo:jhi)) ;swvdr=ZERO ! sw down, visible, direct (W/m^2) allocate(swvdf(ilo:ihi,jlo:jhi)) ;swvdf=ZERO ! sw down, visible, diffuse (W/m^2) allocate(swidr(ilo:ihi,jlo:jhi)) ;swidr=ZERO ! sw down, near IR, direct (W/m^2) allocate(swidf(ilo:ihi,jlo:jhi)) ;swidf=ZERO ! sw down, near IR, diffuse (W/m^2) allocate(flw (ilo:ihi,jlo:jhi)) ;flw =ZERO ! incoming longwave radiation (W/m^2) allocate(frain(ilo:ihi,jlo:jhi)) ;frain=ZERO ! rainfall rate (kg/m^2 s) allocate(fsnow(ilo:ihi,jlo:jhi)) ;fsnow=ZERO ! snowfall rate (kg/m^2 s) ! in from ocean allocate(frzmlt(ilo:ihi,jlo:jhi)) ;frzmlt =ZERO ! freezing/melting potential (W/m^2) allocate(sss (ilo:ihi,jlo:jhi)) ;sss =ZERO ! sea surface salinity (ppt) allocate(sst (ilo:ihi,jlo:jhi)) ;sst =ZERO ! sea surface temperature (C) allocate( Tf (ilo:ihi,jlo:jhi)) ; Tf =ZERO ! freezing temperature (C) allocate( hmix (ilo:ihi,jlo:jhi)) ; hmix =ZERO ! mixed layer depth (m) allocate( qdp (ilo:ihi,jlo:jhi)) ; qdp =ZERO ! deep ocean heat flux (W/m^2) ! out to atmosphere ! note albedos are in ice_albedo.F, Tsfc in ice_state.F ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: & allocate( fsens (ilo:ihi,jlo:jhi)) ;fsens =ZERO ! sensible heat flux (W/m^2) allocate( flat (ilo:ihi,jlo:jhi)) ;flat =ZERO ! latent heat flux (W/m^2) allocate( fswabs(ilo:ihi,jlo:jhi)) ;fswabs =ZERO ! shortwave flux absorbed in ice and ocean (W/m^2) allocate( flwout(ilo:ihi,jlo:jhi)) ;flwout =ZERO ! outgoing longwave radiation (W/m^2) allocate( evap (ilo:ihi,jlo:jhi)) ;evap =ZERO ! evaporative water flux (kg/m^2/s) allocate( Tref (ilo:ihi,jlo:jhi)) ;Tref =ZERO ! 2m atm reference temperature (K) allocate( Qref (ilo:ihi,jlo:jhi)) ;Qref =ZERO ! 2m atm reference sp humidity (kg/kg) ! out to ocean allocate( fresh (ilo:ihi,jlo:jhi)) ;fresh =ZERO ! fresh water flux to ocean (kg/m^2/s) allocate( fsalt (ilo:ihi,jlo:jhi)) ;fsalt =ZERO ! salt flux to ocean (kg/m^2/s) allocate( fhnet (ilo:ihi,jlo:jhi)) ;fhnet =ZERO ! net heat flux to ocean (W/m^2) allocate( fswthru (ilo:ihi,jlo:jhi)) ;fswthru =ZERO ! shortwave penetrating to ocean (W/m^2) ! afm 20150513 for outputs of heat fluxes. ! allocate( fsh(ilo:ihi,jlo:jhi)) ;fsh =ZERO ! sensible heat flux at water surface (W/m^2) ! allocate( flh(ilo:ihi,jlo:jhi)) ;flh =ZERO ! latent heat flux at water surface (W/m^2) ! allocate( flwup(ilo:ihi,jlo:jhi)) ;flwup =ZERO ! outgoing longwave radiation at water surface (W/m^2) ! allocate( flwtot(ilo:ihi,jlo:jhi)) ;flwtot =ZERO ! net longwave radiation at water surface (W/m^2) ! afm 20190807 - bounds fix allocate( fsh(0:MT,jlo:jhi)) ;fsh =ZERO ! sensible heat flux at water surface (W/m^2) allocate( flh(0:MT,jlo:jhi)) ;flh =ZERO ! latent heat flux at water surface (W/m^2) allocate( flwup(0:MT,jlo:jhi)) ;flwup =ZERO ! outgoing longwave radiation at water surface (W/m^2) allocate( flwtot(0:MT,jlo:jhi)) ;flwtot =ZERO ! net longwave radiation at water surface (W/m^2) ! afm 20170125 - hflx ! allocate( flwitot(ilo:ihi,jlo:jhi)) ;flwitot =ZERO ! net longwave radiation at ice surface (W/m^2) ! afm 20190807 - bounds fix allocate( flwitot(0:MT,jlo:jhi)) ;flwitot =ZERO ! net longwave radiation at ice surface (W/m^2) ! diagnostic ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: & allocate(congel(ilo:ihi,jlo:jhi)) ;congel =ZERO! basal ice growth (m/step-->cm/day) allocate(frazil(ilo:ihi,jlo:jhi)) ;frazil =ZERO! frazil ice growth (m/step-->cm/day) allocate(snoice(ilo:ihi,jlo:jhi)) ;snoice =ZERO! snow-ice formation (m/step-->cm/day) allocate(meltt (ilo:ihi,jlo:jhi)) ;meltt =ZERO! top ice melt (m/step-->cm/day) allocate(meltb (ilo:ihi,jlo:jhi)) ;meltb =ZERO! basal ice melt (m/step-->cm/day) allocate(meltl (ilo:ihi,jlo:jhi)) ;meltl =ZERO! lateral ice melt (m/step-->cm/day) allocate(daidtt(ilo:ihi,jlo:jhi)) ;daidtt =ZERO! ice area tendency thermo. (s^-1) allocate(dvidtt(ilo:ihi,jlo:jhi)) ;dvidtt =ZERO! ice volume tendency thermo. (m/s) allocate(mlt_onset(ilo:ihi,jlo:jhi)) ;mlt_onset =ZERO! day of year that sfc melting begins allocate(frz_onset(ilo:ihi,jlo:jhi)) ;frz_onset =ZERO! day of year that freezing begins (congel or frazil) ! NOTE: The following ocean diagnostic fluxes measure ! the same thing as their coupler counterparts but over ! different time intervals. The coupler variables are ! computed from one to_coupler call to the next, whereas ! the diagnostic variables are computed over a time step. allocate(fresh_hist (ilo:ihi,jlo:jhi)) ;fresh_hist =ZERO! fresh water flux to ocean (kg/m^2/s) allocate(fsalt_hist (ilo:ihi,jlo:jhi)) ;fsalt_hist =ZERO! salt flux to ocean (kg/m^2/s) allocate(fhnet_hist (ilo:ihi,jlo:jhi)) ;fhnet_hist =ZERO! net heat flux to ocean (W/m^2) allocate(fswthru_hist(ilo:ihi,jlo:jhi)) ;fswthru_hist=ZERO! shortwave penetrating to ocean (W/m^2) ! internal ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: & allocate(fsw (ilo:ihi,jlo:jhi)) ;fsw =ZERO! ! incoming shortwave radiation (W/m^2) allocate(wind (ilo:ihi,jlo:jhi)) ;wind =ZERO ! wind speed (m/s) allocate(shcoef (ilo:ihi,jlo:jhi)) ;shcoef=ZERO ! transfer coefficient for sensible heat allocate(lhcoef (ilo:ihi,jlo:jhi)) ;lhcoef=ZERO ! transfer coefficient for latent heat !c======================================================================= ! subroutine ice_flux_in_def allocate(cldf(ilo:ihi,jlo:jhi)) ! real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,2) :: & !allocate(fsw_data (ilo:ihi,jlo:jhi,2)) ;fsw_data =ZERO! field values at 2 temporal data points !allocate(cldf_data (ilo:ihi,jlo:jhi,2)) ;cldf_data=ZERO !allocate(fsnow_data(ilo:ihi,jlo:jhi,2)) ;fsnow_data=ZERO !allocate(Tair_data (ilo:ihi,jlo:jhi,2)) ;Tair_data=ZERO !allocate(uatm_data (ilo:ihi,jlo:jhi,2)) ;uatm_data=ZERO !allocate(vatm_data (ilo:ihi,jlo:jhi,2)) ;vatm_data=ZERO !allocate( Qa_data (ilo:ihi,jlo:jhi,2)) ; Qa_data =ZERO !allocate(rhoa_data (ilo:ihi,jlo:jhi,2)) ;rhoa_data=ZERO !allocate(potT_data (ilo:ihi,jlo:jhi,2)) ;potT_data=ZERO !allocate(zlvl_data (ilo:ihi,jlo:jhi,2)) ;zlvl_data=ZERO !allocate( flw_data (ilo:ihi,jlo:jhi,2)) ; flw_data=ZERO !allocate(sst_data (ilo:ihi,jlo:jhi,2)) ;sst_data =ZERO !allocate(sss_data (ilo:ihi,jlo:jhi,2)) ;sss_data =ZERO !======================================================================= ! subroutine ice_mechred_def ! real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi) :: & allocate(dardg1dt(ilo:ihi,jlo:jhi)) ;dardg1dt=ZERO ! rate of fractional area loss by ridging ice (1/s) allocate(dardg2dt(ilo:ihi,jlo:jhi)) ;dardg2dt=ZERO ! rate of fractional area gain by new ridges (1/s) allocate(dvirdgdt(ilo:ihi,jlo:jhi)) ;dvirdgdt=ZERO ! rate of ice volume ridged (m/s) allocate(opening (ilo:ihi,jlo:jhi)) ;opening =ZERO ! rate of opening due to divergence/shear (1/s) ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: & allocate(asum (ilo:ihi,jlo:jhi)) ;asum =ZERO ! sum of total ice and open water area allocate(aksum(ilo:ihi,jlo:jhi)) ;aksum=ZERO ! ratio of area removed to area ridged ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,0:ncat) :: & allocate(apartic(ilo:ihi,jlo:jhi,0:ncat)) ; apartic =ZERO ! participation function; fraction of ridging/ ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,ncat) :: & allocate(hrmin (ilo:ihi,jlo:jhi,ncat)) ;hrmin=ZERO ! minimum ridge thickness allocate(hrmax (ilo:ihi,jlo:jhi,ncat)) ;hrmax=ZERO ! maximum ridge thickness allocate(krdg (ilo:ihi,jlo:jhi,ncat)) ;krdg =ZERO ! mean ridge thickness/thickness of ridging ice allocate(hrexp (ilo:ihi,jlo:jhi,ncat)) ;hrexp=ZERO ! ridge e-folding thickness (krdg_redist = 1) !======================================================================= ! subroutine ice_therm_vertical_def ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,ncat) :: & allocate(hicen_old(ilo:ihi,jlo:jhi,ncat)) ;hicen_old =ZERO ! old ice thickness (m), used by ice_therm_itd ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi) :: & allocate(rside(ilo:ihi,jlo:jhi)) ;rside =ZERO ! fraction of ice that melts laterally !======================================================================= !subroutine ice_therm_itd ! real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi,ncat) :: & allocate(hicen(imt_local,jmt_local,ncat)) ;hicen=ZERO ! ice thickness (m) !======================================================================= ! subroutine ice_grid_def !allocate(TLAT_G(imt_global,jmt_global)) ;TLAT_G=ZERO! latitude of cell center T grid !allocate(TLON_G(imt_global,jmt_global)) ;TLON_G=ZERO! longitude of cell center T grid ! real (kind=dbl_kind), dimension (imt_local,jmt_local) :: & !allocate(dxt(imt_local,jmt_local)) ;dxt=ZERO! width of T-cell through the middle (m) !allocate(dyt(imt_local,jmt_local)) ;dyt=ZERO! height of T-cell through the middle (m) !allocate(dxu(imt_local,jmt_local)) ;dxu=ZERO! width of U-cell through the middle (m) !allocate(dyu(imt_local,jmt_local)) ;dyu=ZERO! height of U-cell through the middle (m) !allocate(HTE(imt_local,jmt_local)) ;HTE=ZERO! length of eastern edge of T-cell (m) !allocate(HTN(imt_local,jmt_local)) ;HTN=ZERO! length of northern edge of T-cell (m) !allocate(HTS(imt_local,jmt_local)) ;HTS=ZERO! length of southern edge of T-cell !allocate(HTW(imt_local,jmt_local)) ;HTW=ZERO! length of western edge of T-cell !allocate(tarea(imt_local,jmt_local)) ;tarea=ZERO! area of T-cell (m^2) !allocate(uarea(imt_local,jmt_local)) ;uarea=ZERO! area of U-cell (m^2) allocate(ULON (imt_local,jmt_local)) ;ULON =ZERO! longitude of velocity pts (radians) allocate(ULAT (imt_local,jmt_local)) ;ULAT =ZERO! latitude of velocity pts (radians) allocate(TLON (imt_local,jmt_local)) ;TLON =ZERO! longitude of temp pts (radians) allocate(TLAT (imt_local,jmt_local)) ;TLAT =ZERO! latitude of temp pts (radians) !allocate(dxhy (imt_local,jmt_local)) ;dxhy =ZERO! 0.5*(HTE - HTE) !allocate(dyhx (imt_local,jmt_local)) ;dyhx =ZERO! 0.5*(HTN - HTN) ! real (kind=dbl_kind), dimension (ilo:;ihi,jlo:jhi) :: & !allocate(cyp (ilo:ihi,jlo:jhi)) ;cyp =ZERO! 1.5*HTE - 0.5*HTE !allocate(cxp (ilo:ihi,jlo:jhi)) ;cxp =ZERO! 1.5*HTN - 0.5*HTN !allocate(cym (ilo:ihi,jlo:jhi)) ;cym =ZERO! 0.5*HTE - 1.5*HTE !allocate(cxm (ilo:ihi,jlo:jhi)) ;cxm =ZERO! 0.5*HTN - 1.5*HTN !allocate(dxt2(ilo:ihi,jlo:jhi)) ;dxt2=ZERO! 0.5*dxt !allocate(dyt2(ilo:ihi,jlo:jhi)) ;dyt2=ZERO! 0.5*dyt !allocate(dxt4(ilo:ihi,jlo:jhi)) ;dxt4=ZERO! 0.25*dxt !allocate(dyt4(ilo:ihi,jlo:jhi)) ;dyt4=ZERO! 0.25*dyt !allocate(tarear (ilo:ihi,jlo:jhi)) ;tarear =ZERO! 1/tarea !allocate(uarear (ilo:ihi,jlo:jhi)) ;uarear =ZERO! 1/uarea !allocate(tinyarea(ilo:ihi,jlo:jhi)) ;tinyarea=ZERO! puny*tarea !allocate(ANGLE (ilo:ihi,jlo:jhi)) ;ANGLE =ZERO! for conversions between POP grid and lat/lon !allocate(ANGLET (ilo:ihi,jlo:jhi)) ;ANGLET =ZERO! ANGLE converted to T-cells !allocate(tarean (ilo:ihi,jlo:jhi)) ;tarean =ZERO! area of NH T-cells !allocate(tareas (ilo:ihi,jlo:jhi)) ;tareas =ZERO! area of SH T-cells ! Masks ! real (kind=dbl_kind), dimension (imt_local,jmt_local) :: & allocate(hm (imt_local,jmt_local)) ; hm =zero! land/boundary mask, thickness (T-cell) !allocate(uvm (imt_local,jmt_local)) ; uvm =ZERO! land/boundary mask, velocity (U-cell) !allocate(mask_n(imt_local,jmt_local)) ; mask_n =ZERO! northern hemisphere !allocate(mask_s(imt_local,jmt_local)) ; mask_s =ZERO! southern hemisphere ! logical (kind=log_kind) :: & allocate(tmask(imt_local,jmt_local)) ;tmask =.true. ! land/boundary mask, thickness (T-cell) allocate(umask(imt_local,jmt_local)) ;umask =.true. ! land/boundary mask, velocity (U-cell) ! allocate(tmask(imt_local,jmt_local)) ;tmask =c1i ! land/boundary mask, thickness (T-cell) ! allocate(umask(imt_local,jmt_local)) ;umask =c1i ! land/boundary mask, velocity (U-cell) !allocate(icetmask(ilo:ihi,jlo:jhi)) ;icetmask =ZERO! ice extent mask (T-cell) !allocate(iceumask(ilo:ihi,jlo:jhi)) ;iceumask =ZERO! ice extent mask (U-cell) !==================================================================== ! my_task, & ! task id for local process ! master_task ! task id for master process my_task=MYID ! processor id master_task = 1 ! master processer id !-------------------------------------------------------------- !==================================================================== ! timing devices ! call ice_timer_clear(-1) ! initialize all timers ! call ice_timer_start(0) ! begin timing entire run ! parameters, grid, parameterizations, etc call init_constants ! pi, etc call input_data ! namelist variables call init_grid ! grid variables ! if (advection == 'remap') ! & call init_remap ! more grid variables ! call init_calendar ! initialize some calendar stuff ! call init_hist ! initialize output history file ! call init_evp ! define evp dynamics parameters, variables call init_flux ! initialize coupler fluxes call init_thermo_vertical ! initialize vertical thermodynamics call init_mechred ! initialize ridging variables call init_itd ! initialize ice thickness distribution ! set the initial ice field base on the latitude ! if (restart) call restartfile ! start from restart data ! call calendar(time) ! determine the initial date !#ifdef coupled ! call init_cpl ! initialize message passing !#else call init_getflux ! initialize forcing data info !#endif Do I=1,M !T Tair(I,1)=273.15+T_air(I) END DO call init_state ! initialize the ice state dtice = dti IF(DBG_SET(DBG_LOG)) write(ipt,*)' ICE time step set to DTI dtice =' , dtice call ALLOC_UVICE call albedos ! albedo based on initial ice distribution IF(DBG_SET(DBG_SBR)) write(ipt,*)'END ice_init_0' return end subroutine ice_init_0 !======================================================================= ! !IROUTINE: from_coupler - input from coupler to sea ice model ! ! !INTERFACE: ! subroutine from_coupler ! !!DESCRIPTION: ! ! Reads input data from coupler to sea ice model ! IMPLICIT NONE integer :: i,j,k,icn,n2 ! local loop indices integer :: II,i1,I2,j1,j2 ! local loop indices real :: gsum, workx, worky REAL(SP), allocatable, dimension (:,:) :: rbuf ! first step initialized the ice variables zlvl(:,:) = C10I uatm(:,:) = 0.0_SP vatm(:,:) = 0.0_SP potT(:,:) = 0.0_SP Tair(:,:) = 0.0_SP Qa(:,:) = 0.0_SP rhoa(:,:) = 1.3_SP swvdr(:,:) = 0.0_SP swvdf(:,:) = 0.0_SP swidr(:,:) = 0.0_SP swidf(:,:) = 0.0_SP flw(:,:) = 0.0_SP frain(:,:) = 0.0_SP fsnow(:,:) = 0.0_SP sst(:,:) = 0.0_SP sss(:,:) = 0.0_SP !uocn(:,:) = 0.0_SP !vocn(:,:) = 0.0_SP !ss_tltx(:,:) = 0.0_SP !ss_tlty(:,:) = 0.0_SP frzmlt(:,:) = 0.0_SP do j=jlo,jhi do i=ilo,ihi !--- ocn states-- sst (i,j) = T1(I,1) !+Tffresh sss (i,j) = S1(I,1) Tf (i,j) = -depressT*sss(i,j) ! freezing temp (C) Tf (i,j) = max(Tf(i,j),-1.85) ! freezing temp (C) !--- atm states- zlvl (i,j) = c10i Tair (i,j) = T_air(I) +Tffresh ! ^0C--->"K" potT (i,j) = Tair(i,j) ! K Qa (i,j) = QA_AIR(I) rhoa (i,j) = 1.3_SP !rhoair !--- ocn states-- ! compute potential to freeze or melt ice frzmlt(i,j) = (Tf(i,j)-sst(i,j))*cp_ocn*rhow*(D(I))*DZ(I,1)/dtice # if !defined (ICE_FRESHWATER) frzmlt(i,j) = min(max(frzmlt(i,j),-c1000),c1000) # endif if (sst(i,j) <= Tf(i,j)) then sst(i,j) = Tf(i,j) # if defined (ICE_FRESHWATER) ! afm 20151112 & EJA 20160921 ! When supercooling (freezing, i.e. frzmlt>0), ! T1 (temperature variable in hydrodynamic modules) is ! set to Tf in subroutine to_coupler. ! The original method only counts surface layer k=1 do k=2,kbm1 IF(T1(I,K) .LT. Tf(I,j)) THEN frzmlt(i,j) = frzmlt(i,j) + (Tf(i,j)-T1(i,k))*cp_ocn*rhow*D(I)*DZ(I,k)/dtice END IF end do # else T1(I,1)=Tf(i,j) # endif endif !--- atm fluxes-- !!====================================================== # if defined (ICE_FRESHWATER) ! afm 20151001 & EJA 20160921 - For solar. For non-solar, shortwave from forcing data is used. fsw(i,j) = swrad_watts(i) # else fsw (i,j) = DSW_AIR(I) # endif cldf(i,j) = CLOUD(I) !!====================================================== fsnow(i,j)=0.0_SP fsnow(i,j) = (QPREC(I)+QEVAP(I))*RHOFRESH !m/s---> kg/m^2/s frain(i,j) = c0i if (Tair(i,j) >= Tffresh) then frain(i,j) = fsnow(i,j) fsnow(i,j) = c0i endif end do end do !!====================================================== !!============================================= !! Special treatment the river mouth ! !! No new ice add at river mouth !! Caused by the salinity changing and the !! freezing potential increasing !!============================================= IF(RIVER_INFLOW_LOCATION == 'node') THEN IF(NUMQBC > 0) THEN DO J=1,NUMQBC I1=INODEQ(J) IF(frzmlt(I1,1)>0) frzmlt(I1,1)=0.0 END DO END IF ELSE IF(RIVER_INFLOW_LOCATION == 'edge') THEN IF(NUMQBC > 0) THEN DO J=1,NUMQBC J1=N_ICELLQ(J,1) J2=N_ICELLQ(J,2) IF(frzmlt(J1,1)>0) frzmlt(J1,1)=0.0 IF(frzmlt(J2,1)>0) frzmlt(J2,1)=0.0 II=ICELLQ(J) DO I1=1,3 IF(isonb(NV(II,I1))==0 ) then I2=NV(II,I1) END IF END DO IF(ISICEN(I2)==1 ) then !IF(frzmlt(J1,1)>0) frzmlt(J1,1)=frzmlt(I2,1) !IF(frzmlt(J2,1)>0) frzmlt(J2,1)=frzmlt(I2,1) END IF END DO END IF END IF !!====================================================== uatm =0.0_SP vatm =0.0_SP wind =0.0_SP DO I=1,M DO I1=1, NTVE(I) Uatm(I,1) =Uatm(I,1)+UUWIND(NBVE(I,I1)) Vatm(I,1) =Vatm(I,1)+VVWIND(NBVE(I,I1)) wind (i,1) =wind(i,1)+ sqrt(UUWIND(NBVE(I,I1))**2 + VVWIND(NBVE(I,I1))**2) ! wind speed, m/s END DO Uatm(I,1) = Uatm(I,1) /float(NTVE(I)) Vatm(I,1) = Vatm(I,1) /float(NTVE(I)) Wind(I,1) = sqrt(Uatm(I,1) * Uatm(I,1)+Vatm(I,1)*Vatm(I,1)) END DO # if defined (MULTIPROCESSOR) allocate(rbuf(0:MT,3)) rbuf =0.0 rbuf(1:M,1)=Uatm(1:M,1) rbuf(1:M,2)=Vatm(1:M,1) rbuf(1:M,3)=Wind(1:M,1) if (PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,3,MYID,NPROCS,rbuf) IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,rbuf) Uatm(1:M,1)=rbuf(1:M,1) Vatm(1:M,1)=rbuf(1:M,2) Wind(1:M,1)=rbuf(1:M,3) deallocate(rbuf) # endif !-----+----------------------------------------------------------------+ ! some parameter correction ! Qa, Fsw , SW-->SWvdr SWVDF SWidr SWidf , FLW, frain, frain call prepare_forcing !scale input forcing data end subroutine from_coupler !======================================================================= !BOP ! ! !IROUTINE: to_coupler - send data from sea ice model to coupler ! ! !INTERFACE: ! subroutine to_coupler ! ! !DESCRIPTION: ! ! Sea ice model to coupler ! use ice_flux integer(kind=int_kind) :: i,j,k,icn,n2 ! local loop indices real (kind=dbl_kind) :: & gsum, workx, worky & ! tmps for converting grid , Tsrf (ilo:ihi,jlo:jhi) & ! surface temperature , tauxa(ilo:ihi,jlo:jhi) & ! atmo/ice stress , tauya(ilo:ihi,jlo:jhi) & , tauxo(ilo:ihi,jlo:jhi) & ! ice/ocean stress , tauyo(ilo:ihi,jlo:jhi) & , ailohi(ilo:ihi,jlo:jhi) ! fractional ice area logical :: flag !!================================================================= REAL(SP) :: SPRO,SPCP,ROSEA ! afm 20171115 for heat flux breakdown consistency-- begin REAL(SP) :: SPRO1(0:MT,jlo:jhi),ROSEA1(0:MT,jlo:jhi) ! afm 20171115 for heat flux breakdown consistency-- end REAL(SP), allocatable, dimension (:,:) :: rbuf REAL (SP) :: Naice ! !!================================================================= !--------ice_ocean--------------------------------------------------- real (kind=dbl_kind), parameter :: & cphm = cp_ocn*rhow*20 !hmix REAL (SP) :: DUVI,ice_ocnDX, ice_ocnDY ! for ice-ocean drag !--------atmo-ocean--------------------------------------------------- real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi) ::& delt &! potential temperature difference (K) , delq &! specific humidity difference (kg/kg) , dummy1, dummy2, dummy3, dummy4 ! dummy arrays real (kind=dbl_kind) :: & TsfK &! surface temperature (K) ! afm for heat flux outputs ! , fsh &! sensible heat flux (W/m^2) ! , flh &! latent heat flux (W/m^2) , swabs !&! surface absorbed shortwave heat flux (W/m^2) ! , flwup ! long-wave upward heat flux (W/m^2) !!==================================================================== flag =.false. SPCP = 4.2174E3_SP ROSEA = 1.023E3_SP SPRO = SPCP*ROSEA ! afm 20171115 for heat flux breakdown consistency-- begin # if defined (HEATING_CALCULATED) ! SPCP = 4.21734E3_SP ! ROSEA = 1.023E3_SP ! SPRO = SPCP*ROSEA !# endif !# if defined (HEATING_CALCULATED_GL) ! SPCP = 4.186E3_SP ! ROSEA = 1.000E3_SP ! SPRO = SPCP*ROSEA ! afm 20180723 -- HEATING_CALCULATED_GL option is merged to HEATING_CALCULATED IF(.NOT. HEATING_FRESHWATER)THEN SPCP=4.2174E3_SP ROSEA = 1.023E3_SP SPRO = SPCP*ROSEA ELSE SPCP= 4.186E3_SP ROSEA = 1.000E3_SP SPRO = SPCP*ROSEA END IF # endif # if defined (HEATING_SOLAR) SPCP = 4.186E3_SP ROSEA1(:,1) = (-.0057_SP)*T1(:,1)*T1(:,1) + 0.0249_SP*T1(:,1) + 1000.0_SP SPRO1 = SPCP*ROSEA1 #endif ! afm 20171115 for heat flux breakdown consistency-- end ! calculate the heat budget at the surface of ocean !!==================================================================== !!==================================================================== ! if(.false.) then ! afm 20160622 & EJA 20160921 -- begin # if defined (HEATING_CALCULATED) || (HEATING_SOLAR) ! Use net heat flux (wtsurf) from COARE or SOLAR. ! Import latent, sensible heat fluxes and net longwave radiation ! for output purpose. ! afm 20190807 - bounds fix flh(1:M,jlo) = HLAT_WATTS(1:M) fsh(1:M,jlo) = HSENS_WATTS(1:M) flwtot(1:M,jlo) = LWRAD_WATTS(1:M) # else ! Calculate net heat flux using CICE's algorithm. ! Overwrite wtsurf. call atmo_boundary_layer (1, 'ocn', sst, & dummy1, dummy2, dummy3, dummy4, delt, delq) do j = jlo,jhi do i = ilo,ihi ! ocean surface temperature in Kelvin TsfK = sst(i,j) + Tffresh ! shortwave radiative flux # if defined (ICE_FRESHWATER) ! afm 20151112 & EJA 20160921 - albedo included in the solar subroutine or forcing file swabs = fsw(i,j) # else swabs = (c1i - albocn) * fsw(i,j) # endif ! longwave radiative flux !flwup = -stefan_boltzmann * TsfK**4 ! afm for heat flux outputs flwup(i,j) = -stefan_boltzmann * TsfK**4 ! downward latent and sensible heat fluxes !flh = lhcoef(i,j) * delq(i,j) !fsh = shcoef(i,j) * delt(i,j) ! afm for heat flux outputs flh(i,j) = lhcoef(i,j) * delq(i,j) fsh(i,j) = shcoef(i,j) * delt(i,j) # if defined (ICE_FRESHWATER) ! afm 20151112 & EJA 20160921 - Bug fix to include emissivity, assume similar to ice surface. ! WTSURF(I) = -(fsh + flh + flwup + flw(i,j)+swabs ) /SPRO !WTSURF(I) = -( fsh + flh + emissivity*( flwup + flw(i,j) ) + swabs ) / SPRO ! for heat flux outputs WTSURF(I) = -( fsh(i,j) + flh(i,j) + emissivity*( flwup(i,j) + flw(i,j)) + swabs ) / SPRO ! afm 20170125 - hflx ! emissivity included in the output FLWTOT(I,J) = emissivity * ( flwup(i,j) + flw(i,j) ) # else WTSURF(I) = -(fsh(i,j) + flh(i,j) + flwup(i,j) + flw(i,j)+swabs ) /SPRO # endif # if !defined (ICE_FRESHWATER) ! afm 20151112 albedo included in the solar subroutine SWRAD(I) = (C1i-albocn)*SWRAD(I) # endif end do end do # endif ! afm 20160622 -- end !!==================================================================== do j = jlo,jhi do i = ilo,ihi if (aice(i,j) > puny) then ! afm 20171115 for heat flux breakdown consistency-- begin # if defined (HEATING_SOLAR) SPRO = SPRO1(i,1) # endif ! afm 20171115 for heat flux breakdown consistency-- end !!========================================================== ! When the ocean is coupled to a sea-ice model, the temperature ! surface boundary condiction switches from a Neuman to a Dirichlet ! boundary condiction: the surface heat flux is no more specified, ! but diagnosed from the setting of the freezing temperature at the ! first ocean level ! the net heat including the second part, tha WTSURF(I) = WTSURF(I)*(c1i-aice(I,1)) & -fhnet(I,J) /SPRO & -fswthru(i,j)/SPRO ! include wssurf correction due to 'fresh' here for salt water application. ! the short wave two parts ! first: the open region ! second: penetrating through the ice to ocean SWRAD(I) = SWRAD(I)*(c1i-aice(I,1)) -fswthru(i,j)/SPRO # if defined (HEATING_CALCULATED) || defined (HEATING_SOLAR) ! afm 20170125 - hflx ! ------ begin flwitot(I,1) = flw(I,1)*emissivity + flwout(I,1) HSENS_WATTS(I) = fsens(I,1)*aice(I,1) + fsh(I,1)*(c1i-aice(i,jlo)) HLAT_WATTS(I) = flat(I,1)*aice(I,1) + flh(I,1)*(c1i-aice(i,jlo)) LWRAD_WATTS(I) = flwitot(I,1)*aice(I,1) + flwtot(I,1)*(c1i-aice(i,jlo)) ! ------ end # endif end if enddo ! i enddo # if defined (HEATING_CALCULATED) || defined (HEATING_SOLAR) ! afm 20170125 - hflx ! ------ begin ! j # if defined (MULTIPROCESSOR) IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,HSENS_WATTS,HLAT_WATTS,LWRAD_WATTS) IF(PAR)CALL AEXCHANGE(NC,MYID,NPROCS,FLWITOT) # endif ! ------ end # endif DO I=1,M IF(T1(I,1) .LT. Tf(I,1)) THEN T1(I,1) = Tf(I,1) END IF DO K=2,KBM1 IF(T1(I,K) .LT. Tf(I,1)) THEN T1(I,K) = Tf(I,1) END IF END DO END DO # if defined (MULTIPROCESSOR) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,T1) # endif CALL N2E3D(T1,T) !Shift to Elements ! !!====================================================== !!====================================================== # if defined (MULTIPROCESSOR) IF(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,WTSURF,SWRAD) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,WTSURF,SWRAD) # endif ! afm 20171115 for heat flux breakdown consistency-- begin # if defined (HEATING_SOLAR) wtsurf_watts = -1.0_SP*wtsurf*SPRO1(:,1) swrad_watts = -1.0_SP*swrad *SPRO1(:,1) # else wtsurf_watts = -1.0_SP*wtsurf*SPRO swrad_watts = -1.0_SP*swrad *SPRO # endif # if defined (MULTIPROCESSOR) IF(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,WTSURF_WATTS,SWRAD_WATTS) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,WTSURF_WATTS,SWRAD_WATTS) # endif !! afm 20171115 -- end DO I=1,M IF(ISICEN(i)==1) THEN ! afm 20180801 moved QEVAP ice fractioning to the non ICE_FRESHWATER flag branch ! QEVAP(I)=(1-aice(i,1))*QEVAP(I) !+EVAP(I,1)/RHOFRESH # if !defined (ICE_FRESHWATER) QEVAP(I)=(1-aice(i,1))*QEVAP(I) !+EVAP(I,1)/RHOFRESH # if !defined (SEMI_IMPLICIT) QPREC(I)=(1-aice(i,1))*QPREC(I)+(fresh(I,1)-evap(I,1))/rhofresh # else QPREC(I)=(1-aice(i,1))*QPREC(I)-evap(I,1)/rhofresh !!yding moved fresh out QPICE(I)=fresh(I,1)/rhofresh !!yding treat fresh seperatelly # endif # if defined (ICE_EMBEDDING) QFLICE(I) = fresh(I,1)/rhofresh !!((RHO1(I,1)+1.0_SP)*1.0E3_SP) QTHICE(I)= vice(I,1)*(RHOI/((RHO1(I,1)+1.0_SP)*1.0E3_SP)) !!yding ICE THICKNESS IN THE WATER mlt_lat(i) = D(i)/QTHICE(i) !!yding for lateral melting adjustment # endif # else ! QPREC(I)=(1-aice(i,1))*QPREC(I)+(fresh(I,1)-evap(I,1))/rhofresh ! afm 20151112 & EJA 20160921 - Change in QPREC calculation. ! no mass flux due to ice melt/freezing ! this causes large omega (WTS) at surface, resulting in unrealistic water ! temperature. Temperature flux due to ice formation/melting still are taken ! account of in wtsurf. ! Note: fresh & evap are already weighted by aicen, and fresh is a summation ! of freshwater mass flux and evap, hence (fresh-evap) in the eq. above. ! EJA test comment out line below to prevent ice-fractioning of precip ! 20180801 afm incorporated the change in 3.2.5 ! QPREC(I)=(1-aice(i,1))*QPREC(I) ! QEVAP(I)=(1-aice(i,1))*QEVAP(I) !+EVAP(I,1)/RHOFRESH # endif ! evap :: flux of vapor, atmos to ice (kg m-2 s-1) !! this is the part lost from the ocean and ice system ! < 0 => sublimation ! > 0 => condensation ! Fresh :: flux of water, ice to ocean (kg/m^2/s) ! negetive ---> ice freezing ! positive ---> ice melting # if !defined (SEMI_IMPLICIT) QPREC2(I) = QPREC(I) QEVAP2(I) = QEVAP(I) ! QPICE2(I) = QPICE(I) !!yding !# if defined (ICE_EMBEDDING) ! QFLICE2(I)= QFLICE(I) !!yding !# endif # endif END IF END DO # if defined (ICE_EMBEDDING) QTHICE_incr = QTHICE - QTHICERK !!yding QTHICERK = QTHICE !!SAVE ICE THICKNESS FOR CURRENT STEP # endif # if defined (MULTIPROCESSOR) if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QPREC,QEVAP) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QPREC,QEVAP) # if !defined (ICE_FRESHWATER) # if defined (SEMI_IMPLICIT) !!yding if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QPICE) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QPICE) !!yding # endif # if defined (ICE_EMBEDDING) if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QFLICE) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QFLICE) if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QTHICE) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QTHICE) if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QTHICERK) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QTHICERK) if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QTHICE_incr) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QTHICE_incr) if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,mlt_lat) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,mlt_lat) # endif # endif # if !defined (SEMI_IMPLICIT) if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QPREC2,QEVAP2) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QPREC2,QEVAP2) !!yding ! if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QPICE2) ! IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QPICE2) !!yding # if defined (ICE_EMBEDDING) ! if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,QFLICE2) ! IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,QFLICE2) # endif # endif # endif end subroutine to_coupler !============================================================== !============================================================== ! Main driver routine for FVCOM and UG-CICE and FVCOM ! ! subroutine ice_1D_2D # if defined (MULTIPROCESSOR) USE MOD_PAR # endif USE BCS USE MOD_OBCS IMPLICIT NONE integer :: i,ii,j,ij,K,Ni, JJ REAL(SP), allocatable, dimension (:,:) :: rbuf REAL(SP), allocatable, dimension (:) :: Irbuf !-------------------------------------------------------------- ! for advection of ice variables !-------------------------------------------------------------- integer (kind=int_kind), parameter :: & narr = 1 + 5*ncat ! number of state variable arrays integer (kind=int_kind) :: & narrays ! counter for number of state variable arrays real (kind=dbl_kind) :: works(imt_local,narr) !works(imt_local,narr) real (kind=dbl_kind) :: worke(imt_local,ntilay) real (kind=dbl_kind) :: TTmin,TTmax integer :: iout integer :: kice !-------------------------------------------------------------- !----------------------------------------------------------------- ! initializations !----------------------------------------------------------------- IF(dbg_set(dbg_SBR)) WRITE(IPT,*) "START ICE_1d_2d" ALLOCATE(AICETMP(0:MT)); AICETMP = 0.0_SP ALLOCATE(IMASS1(0:NT)) ; IMASS1 = 0.0_SP ALLOCATE(IMASS(0:MT)) ; IMASS = 0.0_SP ISICEC_OLD = ISICEC ! ice in last step ISICEC = 0 ISICEN = 0 IMASS =0.0_SP IMASS1 =0.0_SP DO I=1,M !T IF (AICE(I,1) >p001) then IMASS(I) = (RHOI*VICE(I,1) + RHOS*VSNO(I,1)) ! kg/m^2 IF (IMASS(I) >p01) ISICEN(I) =1 END IF END DO # if defined (MULTIPROCESSOR) IF(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,IMASS) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,IMASS) # endif CALL N2E2D(IMASS,IMASS1) DO I=1,M AICETMP(I) = AICE(I,1) END DO AIU =0.0_SP # if defined (MULTIPROCESSOR) if(PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,AICETMP) IF(PAR) CALL AEXCHANGE(NC,MYID,NPROCS,AICETMP) # endif CALL N2E2D(AICETMP,AIU) DO I=1,N IF(sum(ISICEN(NV(I,1:3))) > 0 .and.IMASS1(I)>p01 & .and. AIU(I)>P001) & ISICEC(i)=1 END DO Deallocate(AICETMP) DEALLOCATE(IMASS1,IMASS) !----------------------------------------------------------------- !----------------------------------------------------------------- call from_coupler ! get updated info from flux couple Call thermo_vertical ! thermodynamic growth rates and fluxes call to_coupler ! collect/send data to flux coupler call thermo_itd ! thermodynamics and associated itd changes !!===========================================================================!! # if defined (ONE_D_MODEL) call zap_small_areas ! remove categories with very small areas if (ncat > 1) call rebin ! make sure thicknesses are in bounds call aggregate ! aggregate state variables over categories call albedos ! compute ice albedos return ! NO advection for spinup run # endif !! DYNAMICS FOR SEA ICE DO KICE=1,NDYN_DT !! sub loop for ice model CALL ICE_RUNGE_KUTTA !!=================================================================== !----------------------------------------------------------------- ! Make sure state variables are not too close to zero. !----------------------------------------------------------------- # if defined (ICE_FRESHWATER) ! afm 20160516 not needed. ! commented out as mass loss seems significant ! call check_state # else call check_state # endif !----------------------------------------------------------------- ! fill work arrays with fields to be advected !----------------------------------------------------------------- ! two arrays are used for performance (balance memory/cache vs ! number of bound calls); one array or more than two may perform ! better depending on the machine used, number of processors, etc. ! --tested on SGI R2000, using 4 pes for the ice model under MPI !----------------------------------------------------------------- !!================================================================ works =0.0_SP worke =0.0_SP do i = 1,imt_local do j = 1,jmt_local works(I,1) = min(max(aice0(i,j),0.0),1.0) enddo enddo narrays = 1 do ni=1,ncat do i = 1,imt_local do j = 1,jmt_local ! this is the hice adjust the convergence to hice change works(I,narrays+1) = aicen(i,j,ni) works(I,narrays+2) = vicen(i,j,ni) works(I,narrays+3) = vsnon(i,j,ni) works(I,narrays+4) = -aicen(i,j,ni)*Tsfcn(i,j,ni) works(I,narrays+5) = -(esnon(i,j,ni) & +rhos*Lfresh*vsnon(i,j,ni)) ! for Tsnow<0 enddo enddo ! ! if there were 3 arrays in loop, use 3 instead and change narr narrays = narrays + 5 enddo !----------------------------------------------------------------------------- do k=1,ntilay do i = 1,imt_local do j = 1,jmt_local worke(I,k) = -eicen(i,j,k) enddo enddo enddo if (narr /= narrays) & & write(ipt,*) "Wrong number of arrays in transport bound call" !!-------------------------------------------------------------------------- allocate(Irbuf(0:MT)) !!-------------------------------------------------------------------------- Irbuf=0.0_SP ! afm 20210226 initialize xflux_ice_obc (bug) xflux_ice_obc = 0.0_SP DO K=1,narr Irbuf(1:M)=works(1:M,K) # if defined (MULTIPROCESSOR) if (PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,Irbuf) IF (PAR) CALL AEXCHANGE(NC,MYID,NPROCS,Irbuf) #endif call ICE_ADV_fvcom(Irbuf,UICE2,VICE2) ! - begin afm 20171113 for ice open boundary ------------- ! if vice, save XFLUX to XFLUX_ICE for the open boundary treatment. ! xflux_ice includes contributions from all the categories. IF ( MOD( K, 5 ) == 3 ) then DO I=1,IOBCN XFLUX_ICE_OBC(I) = XFLUX_ICE_OBC(I) + XFLUX_ICE(I_OBC_N(I)) END DO END IF ! - end afm 20171113 for ice open boundary ------------- # if defined (MULTIPROCESSOR) if (PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,Irbuf) IF (PAR) CALL AEXCHANGE(NC,MYID,NPROCS,Irbuf) #endif works(1:M,K)=Irbuf(1:M) END DO !!-------------------------------------------------------------------------- Irbuf=0.0_SP do k=1,ntilay Irbuf(1:M)=worke(1:M,K) # if defined (MULTIPROCESSOR) IF (PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,Irbuf) IF (PAR) CALL AEXCHANGE(NC,MYID,NPROCS,Irbuf) #endif call ICE_ADV_fvcom(Irbuf,UICE2,VICE2) # if defined (MULTIPROCESSOR) IF (PAR) CALL NODE_MATCH(1,NBN,BN_MLT,BN_LOC,BNC,MT,1,MYID,NPROCS,Irbuf) IF (PAR) CALL AEXCHANGE(NC,MYID,NPROCS,Irbuf) #endif worke(1:M,K)=Irbuf(1:M) end do deallocate(Irbuf) !!========================================================================== !----------------------------------------------------------------- ! retrieve advected fields from work array !----------------------------------------------------------------- do i = 1,imt_local do j = 1,jmt_local aice0(i,j) = min(max(works(I,1),C0i),c1i) enddo enddo narrays = 1 ! aice0 is first array do ni=1,ncat do i = 1,imt_local do j = 1,jmt_local aicen(i,j,ni) = min(max(works(I,narrays+1),C0i),c1i) vicen(i,j,ni) = works(I,narrays+2) if ((aicen(i,j,ni) puny) then Tsfcn(i,j,ni) = -works(i,narrays+4)/aicen(i,j,ni) else Tsfcn(i,j,ni) = Tf(i,j) endif if (aicen(i,j,ni) > puny) & esnon(i,j,ni) = -works(i,narrays+5) & & -rhos*Lfresh*vsnon(i,j,ni) ! for Tsnow<0 enddo enddo narrays = narrays + 5 enddo J=1 do i = 1,imt_local do ni=1,ncat # if defined (ICE_FRESHWATER) Tsfcn(i,j,ni) = max( min( Tf(i,j), Tsfcn(i,j,ni)), T_air(i) - 10. ) # endif !IF(Tsfcn(i,j,ni) .le.T_air(I)) Tsfcn(i,j,ni) =T_air(i) enddo enddo !============================================================================= do k=1,ntilay do i=1,imt_local do j=1,jmt_local eicen(i,j,k) = -worke(i,k) enddo enddo enddo 150 continue call check_state call aggregate ! aggregate state variables over categories J=1 do i = 1,imt_local IF (aice(I,j) > 1.0)then do ni=1,ncat aicen(i,j,ni) = aicen(i,j,ni)/aice(i,j) vicen(i,j,ni) = vicen(i,j,ni)/aice(i,j) vsnon(i,j,ni) = vsnon(i,j,ni)/aice(i,j) esnon(i,j,ni) = esnon(i,j,ni)/aice(i,j) do k = 1, nilyr eicen(i,j,ilyr1(ni)+k-1) =eicen(i,j,ilyr1(ni)+k-1)/aice(i,j) end do enddo aice0 (i,j)=0.0 aice(i,j)=1.0 END IF enddo ! !! begin afm 20210225 for ice river outflow j=1 IF(RIVER_INFLOW_LOCATION == 'node') THEN IF(NUMQBC > 0) THEN DO I=1,NUMQBC IF(qdis(I) < 0.0_SP) THEN ! IF THE FLOW IS OUT OF THE DOMAIN ! neighboring node is not defined for river. ! just remove when outflowing aice(INODEQ(I),j) = 0.0_SP vice(INODEQ(I),j) = 0.0_SP vicen(INODEQ(I),j,:) = 0.0_DP aicen(INODEQ(I),j,:) = 0.0_DP vsnon(INODEQ(I),j,:) = 0.0_SP esnon(INODEQ(I),j,:) = 0.0_SP eicen(INODEQ(I),j,:) = 0.0_SP aice0(INODEQ(I),j) = 1.0_SP ENDIF ENDDO ENDIF ELSE IF(RIVER_INFLOW_LOCATION == 'edge') THEN IF(NUMQBC > 0) THEN DO I=1,NUMQBC IF(qdis(I) < 0.0_SP) THEN ! IF THE FLOW IS OUT OF THE DOMAIN ! neighboring node is not defined for river. ! just remove when outflowing do k = 1,2 aice(N_ICELLQ(I,k),j) = 0.0_SP vice(N_ICELLQ(I,k),j) = 0.0_SP vicen(N_ICELLQ(I,k),j,:) = 0.0_DP aicen(N_ICELLQ(I,k),j,:) = 0.0_DP vsnon(N_ICELLQ(I,k),j,:) = 0.0_SP esnon(N_ICELLQ(I,k),j,:) = 0.0_SP eicen(N_ICELLQ(I,k),j,:) = 0.0_SP aice0(N_ICELLQ(I,k),j) = 1.0_SP end do ENDIF ENDDO ENDIF ENDIF ! end afm 20210225 for ice river outflow ! begin afm 20210414 limit vice by water depth # if defined (ICE_FRESHWATER) J=1 do i = 1,imt_local IF (VICE(I,J) >= D(I) ) then print*, 'vice > water depth, scale vice. d=h+el, local node:', D(I),VICE(I,J),i do ni=1,ncat vicen(i,j,ni) = vicen(i,j,ni)*D(I)/vice(i,j) aicen(i,j,ni) = aicen(i,j,ni)*D(I)/vice(i,j) vsnon(i,j,ni) = vsnon(i,j,ni)*D(I)/vice(i,j) esnon(i,j,ni) = esnon(i,j,ni)*D(I)/vice(i,j) do k = 1, nilyr eicen(i,j,ilyr1(ni)+k-1)=eicen(i,j,ilyr1(ni)+k-1)*D(I)/vice(i,j) end do enddo ! comment out. aggregation will be done later. ! aice0 (i,j)=0.0 ! aice(i,j)=1.0 END IF enddo # endif ! end afm 20210414 limit vice by water depth call ridge_ice(Delta,divu) ! ridging ENDDO !! KTICE sub loop for ice dynamics call zap_small_areas ! remove categories with very small areas if (ncat > 1) call rebin ! make sure thicknesses are in bounds call aggregate ! aggregate state variables over categories call albedos ! compute ice albedos ! hhu ice bc 2/8/2017 call bcond_gcn(4) ! hhu IF(dbg_set(dbg_SBR)) WRITE(IPT,*) "END ICE_1d_2d" return end subroutine ice_1D_2D !=======================================================================! # endif END MODULE MOD_ICE