318 lines
16 KiB
Fortran
Executable File
318 lines
16 KiB
Fortran
Executable File
MODULE icestp
|
|
!!======================================================================
|
|
!! *** MODULE icestp ***
|
|
!! sea ice : Master routine for all the sea ice model
|
|
!!=====================================================================
|
|
!!
|
|
!! The sea ice model SI3 (Sea Ice modelling Integrated Initiative),
|
|
!! aka Sea Ice cube for its nickname
|
|
!!
|
|
!! is originally based on LIM3, developed in Louvain-la-Neuve by:
|
|
!! * Martin Vancoppenolle (UCL-ASTR, Belgium)
|
|
!! * Sylvain Bouillon (UCL-ASTR, Belgium)
|
|
!! * Miguel Angel Morales Maqueda (NOC-L, UK)
|
|
!! thanks to valuable earlier work by
|
|
!! * Thierry Fichefet
|
|
!! * Hugues Goosse
|
|
!! thanks also to the following persons who contributed
|
|
!! * Gurvan Madec, Claude Talandier, Christian Ethe (LOCEAN, France)
|
|
!! * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany)
|
|
!! * Bill Lipscomb (LANL), Cecilia Bitz (UWa) and Elisabeth Hunke (LANL), USA.
|
|
!!
|
|
!! SI3 has been made possible by a handful of persons who met as working group
|
|
!! (from France, Belgium, UK and Italy)
|
|
!! * Clement Rousset, Martin Vancoppenolle & Gurvan Madec (LOCEAN, France)
|
|
!! * Matthieu Chevalier & David Salas (Meteo France, France)
|
|
!! * Gilles Garric (Mercator Ocean, France)
|
|
!! * Thierry Fichefet & Francois Massonnet (UCL, Belgium)
|
|
!! * Ed Blockley & Jeff Ridley (Met Office, UK)
|
|
!! * Danny Feltham & David Schroeder (CPOM, UK)
|
|
!! * Yevgeny Aksenov (NOC, UK)
|
|
!! * Paul Holland (BAS, UK)
|
|
!! * Dorotea Iovino (CMCC, Italy)
|
|
!!======================================================================
|
|
!! History : 4.0 ! 2018 (C. Rousset) Original code SI3
|
|
!!----------------------------------------------------------------------
|
|
#if defined key_si3
|
|
!!----------------------------------------------------------------------
|
|
!! 'key_si3' SI3 sea-ice model
|
|
!!----------------------------------------------------------------------
|
|
!! ice_stp : sea-ice model time-stepping and update ocean SBC over ice-covered area
|
|
!! ice_init : initialize sea-ice
|
|
!!----------------------------------------------------------------------
|
|
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 phycst ! Define parameters for the routines
|
|
USE eosbn2 ! equation of state
|
|
USE sbc_oce ! Surface boundary condition: ocean fields
|
|
USE sbc_ice ! Surface boundary condition: ice fields
|
|
!
|
|
USE icesbc ! sea-ice: Surface boundary conditions
|
|
USE iceupdate ! sea-ice: sea surface boundary condition update
|
|
USE icewri ! sea-ice: outputs
|
|
!
|
|
!
|
|
USE in_out_manager ! I/O manager
|
|
USE iom ! I/O manager library
|
|
USE lib_mpp ! MPP library
|
|
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
|
|
USE timing ! Timing
|
|
USE prtctl ! Print control
|
|
|
|
IMPLICIT NONE
|
|
PRIVATE
|
|
|
|
PUBLIC ice_stp ! called by sbcmod.F90
|
|
PUBLIC ice_init ! called by sbcmod.F90
|
|
|
|
!! * Substitutions
|
|
# include "do_loop_substitute.h90"
|
|
!!----------------------------------------------------------------------
|
|
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
|
|
!! $Id: icestp.F90 13655 2020-10-21 14:15:13Z laurent $
|
|
!! Software governed by the CeCILL license (see ./LICENSE)
|
|
!!----------------------------------------------------------------------
|
|
CONTAINS
|
|
|
|
SUBROUTINE ice_stp( kt, Kbb, Kmm, 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) :: Kbb, Kmm ! ocean time level indices
|
|
INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk, or Pure Coupled)
|
|
!
|
|
INTEGER :: jl ! dummy loop index
|
|
!!----------------------------------------------------------------------
|
|
!
|
|
IF( ln_timing ) CALL timing_start('ice_stp')
|
|
!
|
|
! !-----------------------!
|
|
IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! --- Ice time step --- !
|
|
! !-----------------------!
|
|
!
|
|
kt_ice = kt ! -- Ice model time step
|
|
!
|
|
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) )
|
|
!
|
|
!
|
|
!------------------------------------------------!
|
|
! --- 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_sbc_tau( kt, ksbc, utau_ice, vtau_ice )
|
|
!
|
|
!------------------------------------------------------!
|
|
! --- Thermodynamical coupling with the atmosphere --- !
|
|
!------------------------------------------------------!
|
|
! It provides the following fields used in sea ice model:
|
|
! 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_sbc_flx( kt, ksbc )
|
|
!----------------------------!
|
|
! --- ice thermodynamics --- !
|
|
!----------------------------!
|
|
|
|
fr_i (:,:) = at_i(:,:) !#LB...
|
|
|
|
!!#LB: lines stolen from "ice_update_flx()@iceupdate.F90"
|
|
!!=============================================================
|
|
|
|
! --- case we bypass ice thermodynamics --- !
|
|
!IF( .NOT. ln_icethd ) THEN ! we suppose ice is impermeable => ocean is isolated from atmosphere
|
|
qt_atm_oi (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
|
|
qt_oce_ai (:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)
|
|
emp_ice (:,:) = 0._wp
|
|
qemp_ice (:,:) = 0._wp
|
|
qevap_ice (:,:,:) = 0._wp
|
|
!ENDIF
|
|
|
|
! output all fluxes
|
|
!------------------
|
|
|
|
! --- mass fluxes [kg/m2/s] --- !
|
|
IF( iom_use('emp_oce' ) ) CALL iom_put( 'emp_oce', emp_oce ) ! emp over ocean (taking into account the snow blown away from the ice)
|
|
IF( iom_use('emp_ice' ) ) CALL iom_put( 'emp_ice', emp_ice ) ! emp over ice (taking into account the snow blown away from the ice)
|
|
! --- heat fluxes [W/m2] --- !
|
|
IF( iom_use('qsr_oce' ) ) CALL iom_put( 'qsr_oce' , qsr_oce * ( 1._wp - at_i_b ) ) ! solar flux at ocean surface
|
|
IF( iom_use('qns_oce' ) ) CALL iom_put( 'qns_oce' , qns_oce * ( 1._wp - at_i_b ) + qemp_oce ) ! non-solar flux at ocean surface
|
|
IF( iom_use('qns_ice' ) ) CALL iom_put( 'qns_ice' , SUM( qns_ice * a_i_b, dim=3 ) + qemp_ice ) ! non-solar flux at ice surface
|
|
IF( iom_use('qsr_ice' ) ) CALL iom_put( 'qsr_ice' , SUM( qsr_ice * a_i_b, dim=3 ) )
|
|
IF( iom_use('qt_oce' ) ) CALL iom_put( 'qt_oce' , ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce )
|
|
IF( iom_use('qt_ice' ) ) CALL iom_put( 'qt_ice' , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 ) + qemp_ice )
|
|
IF( iom_use('qemp_oce' ) ) CALL iom_put( 'qemp_oce' , qemp_oce ) ! Downward Heat Flux from E-P over ocean
|
|
IF( iom_use('qemp_ice' ) ) CALL iom_put( 'qemp_ice' , qemp_ice ) ! Downward Heat Flux from E-P over ice
|
|
|
|
|
|
!!=============================================================
|
|
|
|
!
|
|
CALL ice_wri( kt ) ! -- Ice outputs
|
|
!
|
|
!IF( lrst_ice ) CALL ice_rst_write( kt ) ! -- Ice restart file
|
|
!
|
|
!IF( ln_icectl ) CALL ice_ctl( kt ) ! -- Control checks
|
|
!
|
|
ENDIF ! End sea-ice time step only
|
|
|
|
!-------------------------!
|
|
! --- Ocean time step --- !
|
|
!-------------------------!
|
|
!
|
|
IF( ln_timing ) CALL timing_stop('ice_stp')
|
|
!
|
|
END SUBROUTINE ice_stp
|
|
|
|
|
|
SUBROUTINE ice_init( Kbb, Kmm, Kaa )
|
|
!!----------------------------------------------------------------------
|
|
!! *** ROUTINE ice_init ***
|
|
!!
|
|
!! ** purpose : Initialize sea-ice parameters
|
|
!!----------------------------------------------------------------------
|
|
INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
|
|
!
|
|
INTEGER :: ierr
|
|
!!----------------------------------------------------------------------
|
|
IF(lwp) WRITE(numout,*)
|
|
IF(lwp) WRITE(numout,*) 'Sea Ice Model: SI3 (Sea Ice modelling Integrated Initiative)'
|
|
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
|
|
IF(lwp) WRITE(numout,*)
|
|
IF(lwp) WRITE(numout,*) 'ice_init: Arrays allocation & Initialization of all routines & init state'
|
|
IF(lwp) WRITE(numout,*) '~~~~~~~~'
|
|
!
|
|
! ! Load the reference and configuration namelist files and open namelist output file
|
|
CALL load_nml( numnam_ice_ref, 'namelist_ice_ref', numout, lwm )
|
|
CALL load_nml( numnam_ice_cfg, 'namelist_ice_cfg', numout, lwm )
|
|
IF(lwm) CALL ctl_opn( numoni , 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
|
|
!
|
|
CALL par_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 boundary conditions
|
|
ierr = ierr + ice1D_alloc () ! thermodynamics
|
|
!
|
|
CALL mpp_sum( 'icestp', ierr )
|
|
IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
|
|
!
|
|
!
|
|
CALL ice_sbc_init ! set ice-ocean and ice-atm. coupling parameters
|
|
!
|
|
fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction
|
|
!
|
|
IF( ln_rstart ) CALL iom_close( numrir ) ! close input ice restart file
|
|
!
|
|
END SUBROUTINE ice_init
|
|
|
|
|
|
SUBROUTINE par_init
|
|
!!-------------------------------------------------------------------
|
|
!! *** ROUTINE par_init ***
|
|
!!
|
|
!! ** Purpose : Definition generic parameters for ice model
|
|
!!
|
|
!! ** Method : Read namelist and check the parameter
|
|
!! values called at the first timestep (nit000)
|
|
!!
|
|
!! ** input : Namelist nampar
|
|
!!-------------------------------------------------------------------
|
|
INTEGER :: ios ! Local integer
|
|
!!
|
|
NAMELIST/nampar/ jpl, nlay_i, nlay_s, ln_virtual_itd, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s, &
|
|
& cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir
|
|
!!-------------------------------------------------------------------
|
|
!
|
|
READ ( numnam_ice_ref, nampar, IOSTAT = ios, ERR = 901)
|
|
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampar in reference namelist' )
|
|
READ ( numnam_ice_cfg, nampar, IOSTAT = ios, ERR = 902 )
|
|
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nampar in configuration namelist' )
|
|
IF(lwm) WRITE( numoni, nampar )
|
|
!
|
|
IF(lwp) THEN ! control print
|
|
WRITE(numout,*)
|
|
WRITE(numout,*) ' par_init: ice parameters shared among all the routines'
|
|
WRITE(numout,*) ' ~~~~~~~~'
|
|
WRITE(numout,*) ' Namelist nampar: '
|
|
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 param for jpl=1 (T) or not (F) ln_virtual_itd = ', ln_virtual_itd
|
|
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
|
|
! !--- change max ice concentration for roundoff errors
|
|
rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 )
|
|
rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 )
|
|
! !--- check consistency
|
|
IF ( jpl > 1 .AND. ln_virtual_itd ) THEN
|
|
ln_virtual_itd = .FALSE.
|
|
IF(lwp) WRITE(numout,*)
|
|
IF(lwp) WRITE(numout,*) ' ln_virtual_itd forced to false as jpl>1, no need with multiple categories to emulate them'
|
|
ENDIF
|
|
!
|
|
IF( ln_cpl .AND. nn_cats_cpl /= 1 .AND. nn_cats_cpl /= jpl ) THEN
|
|
CALL ctl_stop( 'STOP', 'par_init: in coupled mode, nn_cats_cpl should be either 1 or jpl' )
|
|
ENDIF
|
|
!
|
|
rDt_ice = REAL(nn_fsbc) * rn_Dt !--- sea-ice timestep and its inverse
|
|
r1_Dt_ice = 1._wp / rDt_ice
|
|
IF(lwp) WRITE(numout,*)
|
|
IF(lwp) WRITE(numout,*) ' ice timestep rDt_ice = nn_fsbc*rn_Dt = ', 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 par_init
|
|
|
|
#else
|
|
!!----------------------------------------------------------------------
|
|
!! Default option Dummy module NO SI3 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
|