MODULE icestp !!====================================================================== !! *** MODULE icestp *** !! Surface module : update the ocean surface boundary condition over ice !! & covered area using LIM sea-ice model !! Sea-Ice model : LIM-3 Sea ice model time-stepping !!===================================================================== !! History : 2.0 ! 2006-12 (M. Vancoppenolle) Original code !! 3.0 ! 2008-02 (C. Talandier) Surface module from icestp.F90 !! - ! 2008-04 (G. Madec) sltyle and ice_ctl routine !! 3.3 ! 2010-11 (G. Madec) ice-ocean stress always computed at each ocean time-step !! 3.4 ! 2011-01 (A Porter) dynamical allocation !! - ! 2012-10 (C. Rousset) add ice_dia !! 3.6 ! 2014-07 (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface !! 4.0 ! 2016-06 (L. Brodeau) new unified bulk routine (based on AeroBulk) !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM 3.0 sea-ice model !!---------------------------------------------------------------------- !! ice_stp : sea-ice model time-stepping and update ocean surf. boundary cond. over ice-covered area !! ice_init : !! ice_run_init : !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE c1d ! 1D vertical configuration USE ice ! sea-ice: variables USE ice1D ! sea-ice: thermodynamical 1D variables ! USE sbc_oce ! Surface boundary condition: ocean fields USE sbc_ice ! Surface boundary condition: ice fields USE iceforcing ! sea-ice: Surface boundary condition !!gm why not icesbc module name ! USE phycst ! Define parameters for the routines USE eosbn2 ! equation of state USE icerhg ! sea-ice: rheology USE iceadv ! sea-ice: advection USE icedyn ! sea-ice: dynamics USE icethd ! sea-ice: thermodynamics USE icerdgrft ! sea-ice: ridging/rafting USE iceupdate ! sea-ice: sea surface boundary condition update USE icedia ! sea-ice: budget diagnostics USE icewri ! sea-ice: outputs USE icerst ! sea-ice: restarts USE icecor ! sea-ice: corrections USE icevar ! sea-ice: operations USE icectl ! sea-ice: control ! MV MP 2016 USE limmp ! sea-ice: melt ponds ! END MV MP 2016 USE iceistate ! sea-ice: initial state USE icethd_sal ! sea-ice: thermodynamics and salinity USE iceitd ! sea-ice: remapping thickness distribution USE icealb ! sea-ice: albedo ! USE bdy_oce , ONLY : ln_bdy ! flag for bdy USE bdyice ! unstructured open boundary data for sea-ice # if defined key_agrif USE agrif_ice USE agrif_lim3_update USE agrif_lim3_interp # endif ! USE in_out_manager ! I/O manager USE iom ! I/O manager library USE prtctl ! Print control USE lib_fortran ! USE lbclnk ! lateral boundary condition - MPP link USE lib_mpp ! MPP library USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC ice_stp ! called by sbcmod.F90 PUBLIC ice_init ! called by sbcmod.F90 !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2017) !! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_stp( kt, ksbc ) !!--------------------------------------------------------------------- !! *** ROUTINE ice_stp *** !! !! ** Purpose : sea-ice model time-stepping and update ocean surface !! boundary condition over ice-covered area !! !! ** Method : ice model time stepping !! - call the ice dynamics routine !! - call the ice advection/diffusion routine !! - call the ice thermodynamics routine !! - call the routine that computes mass and !! heat fluxes at the ice/ocean interface !! - save the outputs !! - save the outputs for restart when necessary !! !! ** Action : - time evolution of the LIM sea-ice model !! - update all sbc variables below sea-ice: !! utau, vtau, taum, wndm, qns , qsr, emp , sfx !!--------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time step INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk, or Pure Coupled) ! INTEGER :: jl ! dummy loop index !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('ice_stp') ! ! !-----------------------! IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! --- Ice time step --- ! ! !-----------------------! ! kt_ice = kt ! -- Ice model time step ! # if defined key_agrif IF( .NOT. Agrif_Root() ) lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1 # endif u_oce(:,:) = ssu_m(:,:) ! -- mean surface ocean current v_oce(:,:) = ssv_m(:,:) ! CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) ! -- freezing temperature [Kelvin] (set to rt0 over land) t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) ! CALL store_fields ! -- Store now ice values ! !------------------------------------------------! ! --- Dynamical coupling with the atmosphere --- ! !------------------------------------------------! ! It provides the following fields used in sea ice model: ! utau_ice, vtau_ice = surface ice stress [N/m2] !------------------------------------------------! CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice ) !-------------------------------------! ! --- ice dynamics and advection --- ! !-------------------------------------! CALL ice_diag0 ! set diag of mass, heat and salt fluxes to 0 CALL ice_rst_opn( kt ) ! Open Ice restart file (if necessary) ! IF( ln_icedyn .AND. .NOT.lk_c1d ) CALL ice_dyn( kt ) ! -- Ice dynamics ! ! !== lateral boundary conditions ==! #if defined key_agrif IF( .NOT. Agrif_Root() ) CALL agrif_interp_lim3('T') ! -- AGRIF #endif IF( ln_icethd .AND. ln_bdy ) CALL bdy_ice( kt ) ! -- bdy ice thermo ! ! ! !== previous lead fraction and ice volume for flux calculations ! CALL ice_var_glo2eqv ! ht_i and ht_s for ice albedo calculation CALL ice_var_agg(1) ! at_i for coupling CALL store_fields ! Store now ice values !------------------------------------------------------! ! --- Thermodynamical coupling with the atmosphere --- ! !------------------------------------------------------! ! It provides the following fields used in sea ice model: ! fr1_i0 , fr2_i0 = 1sr & 2nd fraction of qsr penetration in ice [%] ! emp_oce , emp_ice = E-P over ocean and sea ice [Kg/m2/s] ! sprecip = solid precipitation [Kg/m2/s] ! evap_ice = sublimation [Kg/m2/s] ! qsr_tot , qns_tot = solar & non solar heat flux (total) [W/m2] ! qsr_ice , qns_ice = solar & non solar heat flux over ice [W/m2] ! dqns_ice = non solar heat sensistivity [W/m2] ! qemp_oce, qemp_ice, = sensible heat (associated with evap & precip) [W/m2] ! qprec_ice, qevap_ice !------------------------------------------------------! CALL ice_forcing_flx( kt, ksbc ) !----------------------------! ! --- ice thermodynamics --- ! !----------------------------! IF( ln_icethd ) CALL ice_thd( kt ) ! -- Ice thermodynamics ! MV MP 2016 IF ( ln_pnd ) CALL lim_mp( kt ) ! -- Melt ponds ! END MV MP 2016 IF( ln_icethd ) CALL ice_cor( kt , 2 ) ! -- Corrections ! --- # if defined key_agrif IF( .NOT. Agrif_Root() ) CALL agrif_update_lim3( kt ) # endif CALL ice_var_glo2eqv ! necessary calls (at least for coupling) CALL ice_var_agg( 2 ) ! necessary calls (at least for coupling) ! !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work !! moreover it should only be called at the update frequency (as in agrif_lim3_update.F90) !!# if defined key_agrif !! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid() !!# endif CALL ice_update_flx( kt ) ! -- Update ocean surface mass, heat and salt fluxes !!# if defined key_agrif !! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid() !!# endif IF( ln_icediahsb ) CALL ice_dia( kt ) ! -- Diagnostics and outputs ! CALL ice_wri( 1 ) ! -- Ice outputs ! IF( kt == nit000 .AND. ln_rstart ) & !!gm I don't understand the ln_rstart, if needed, add a comment, please & CALL iom_close( numrir ) ! close input ice restart file ! IF( lrst_ice ) CALL ice_rst_write( kt ) ! -- Ice restart file ! IF( ln_icectl ) CALL ice_ctl( kt ) ! alerts in case of model crash ! ENDIF ! End sea-ice time step only !-------------------------! ! --- Ocean time step --- ! !-------------------------! ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere ! using before instantaneous surf. currents IF( ln_icedyn ) CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) ) !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! ! IF( nn_timing == 1 ) CALL timing_stop('ice_stp') ! END SUBROUTINE ice_stp SUBROUTINE ice_init !!---------------------------------------------------------------------- !! *** ROUTINE ice_init *** !! !! ** purpose : Allocate all the dynamic arrays of the LIM-3 modules !!---------------------------------------------------------------------- INTEGER :: ji, jj, ierr !!---------------------------------------------------------------------- IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'ice_init: Arrays allocation & Initialization off all routines & init state' IF(lwp) WRITE(numout,*) '~~~~~~~~' ! ! ! Open the reference and configuration namelist files and namelist output file CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 ) ! CALL ice_run_init ! set some ice run parameters ! ! ! Allocate the ice arrays (sbc_ice already allocated in sbc_init) ierr = ice_alloc () ! ice variables ierr = ierr + sbc_ice_alloc () ! surface forcing ierr = ierr + ice1D_alloc () ! thermodynamics ierr = ierr + ice_rdgrft_alloc () ! ridging/rafting ! IF( lk_mpp ) CALL mpp_sum( ierr ) IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays') ! CALL ice_itd_init ! ice thickness distribution initialization ! IF( ln_icedyn ) THEN CALL ice_dyn_init ! set ice dynamics parameters CALL ice_rdgrft_init ! set ice ridging/rafting parameters CALL ice_rhg_init ! set ice rheology parameters CALL ice_adv_init ! set ice advection parameters ENDIF IF( ln_icethd ) THEN CALL ice_thd_init ! set ice thermodynics parameters CALL ice_thd_sal_init ! set ice salinity parameters ENDIF ! MV MP 2016 CALL lim_mp_init ! set melt ponds parameters ! END MV MP 2016 ! ! Initial sea-ice state IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst CALL ice_istate_init CALL ice_istate ELSE ! start from a restart file CALL ice_rst_read ENDIF CALL ice_var_agg(2) CALL ice_var_glo2eqv ! CALL ice_update_init ! ice surface boundary condition ! CALL ice_alb_init ! ice surface albedo ! CALL ice_dia_init ! initialization for diags ! fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu ! ! ! set max concentration in both hemispheres WHERE( gphit(:,:) > 0._wp ) ; rn_amax_2d(:,:) = rn_amax_n ! NH ELSEWHERE ; rn_amax_2d(:,:) = rn_amax_s ! SH END WHERE ! END SUBROUTINE ice_init SUBROUTINE ice_run_init !!------------------------------------------------------------------- !! *** ROUTINE ice_run_init *** !! !! ** Purpose : Definition some run parameter for ice model !! !! ** Method : Read the namice_run namelist and check the parameter !! values called at the first timestep (nit000) !! !! ** input : Namelist namice_run !!------------------------------------------------------------------- INTEGER :: ios ! Local integer output status for namelist read NAMELIST/namice_run/ jpl, nlay_i, nlay_s, nn_monocat, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s, & & cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir !!------------------------------------------------------------------- ! REWIND( numnam_ice_ref ) ! Namelist namice_run in reference namelist : Parameters for ice READ ( numnam_ice_ref, namice_run, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_run in reference namelist', lwp ) REWIND( numnam_ice_cfg ) ! Namelist namice_run in configuration namelist : Parameters for ice READ ( numnam_ice_cfg, namice_run, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_run in configuration namelist', lwp ) IF(lwm) WRITE ( numoni, namice_run ) ! IF(lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice' WRITE(numout,*) ' ~~~~~~' WRITE(numout,*) ' Namelist namice_run : ' WRITE(numout,*) ' number of ice categories jpl = ', jpl WRITE(numout,*) ' number of ice layers nlay_i = ', nlay_i WRITE(numout,*) ' number of snow layers nlay_s = ', nlay_s WRITE(numout,*) ' virtual ITD mono-category param (1-4) or not (0) nn_monocat = ', nn_monocat WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_icedyn = ', ln_icedyn WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_icethd = ', ln_icethd WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s ENDIF ! ! !--- check consistency IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN nn_monocat = 0 IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen' ENDIF IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' ) ENDIF ! IF( ln_bdy .AND. ln_icediachk ) CALL ctl_warn('ice_run_init: online conservation check does not work with BDY') ! rdt_ice = REAL(nn_fsbc) * rdt !--- sea-ice timestep and inverse r1_rdtice = 1._wp / rdt_ice IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice ! r1_nlay_i = 1._wp / REAL( nlay_i, wp ) !--- inverse of nlay_i and nlay_s r1_nlay_s = 1._wp / REAL( nlay_s, wp ) ! END SUBROUTINE ice_run_init SUBROUTINE store_fields !!---------------------------------------------------------------------- !! *** ROUTINE store_fields *** !! !! ** purpose : store ice variables at "before" time step !!---------------------------------------------------------------------- INTEGER :: ji, jj, jl ! dummy loop index !!---------------------------------------------------------------------- ! a_i_b (:,:,:) = a_i (:,:,:) ! ice area v_i_b (:,:,:) = v_i (:,:,:) ! ice volume v_s_b (:,:,:) = v_s (:,:,:) ! snow volume smv_i_b(:,:,:) = smv_i(:,:,:) ! salt content oa_i_b (:,:,:) = oa_i (:,:,:) ! areal age content e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy WHERE( a_i_b(:,:,:) >= epsi20 ) ht_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:) ! ice thickness ht_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:) ! snw thickness ELSEWHERE ht_i_b(:,:,:) = 0._wp ht_s_b(:,:,:) = 0._wp END WHERE ! ice velocities & total concentration at_i_b(:,:) = SUM( a_i_b(:,:,:), dim=3 ) u_ice_b(:,:) = u_ice(:,:) v_ice_b(:,:) = v_ice(:,:) ! END SUBROUTINE store_fields SUBROUTINE ice_diag0 !!---------------------------------------------------------------------- !! *** ROUTINE ice_diag0 *** !! !! ** purpose : set ice-ocean and ice-atm. fluxes to zeros at the beggining !! of the time step !!---------------------------------------------------------------------- INTEGER :: ji, jj ! dummy loop index !!---------------------------------------------------------------------- sfx (:,:) = 0._wp ; sfx_bri(:,:) = 0._wp ; sfx_lam(:,:) = 0._wp sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp ! wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp wfx_spr(:,:) = 0._wp ; wfx_lam(:,:) = 0._wp wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp wfx_snw_sni(:,:) = 0._wp ! MV MP 2016 wfx_pnd(:,:) = 0._wp ! END MV MP 2016 hfx_thd(:,:) = 0._wp ; hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp hfx_err_dif(:,:) = 0._wp wfx_err_sub(:,:) = 0._wp ! afx_tot(:,:) = 0._wp ; ! diag_heat(:,:) = 0._wp ; diag_smvi(:,:) = 0._wp diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp ! SIMIP diagnostics diag_fc_bo(:,:) = 0._wp ; diag_fc_su(:,:) = 0._wp tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here) END SUBROUTINE ice_diag0 #else !!---------------------------------------------------------------------- !! Default option Dummy module NO LIM 3.0 sea-ice model !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_stp ( kt, ksbc ) ! Dummy routine INTEGER, INTENT(in) :: kt, ksbc WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt END SUBROUTINE ice_stp SUBROUTINE ice_init ! Dummy routine END SUBROUTINE ice_init #endif !!====================================================================== END MODULE icestp