[13655] | 1 | MODULE icestp |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE icestp *** |
---|
| 4 | !! sea ice : Master routine for all the sea ice model |
---|
| 5 | !!===================================================================== |
---|
| 6 | !! |
---|
| 7 | !! The sea ice model SI3 (Sea Ice modelling Integrated Initiative), |
---|
| 8 | !! aka Sea Ice cube for its nickname |
---|
| 9 | !! |
---|
| 10 | !! is originally based on LIM3, developed in Louvain-la-Neuve by: |
---|
| 11 | !! * Martin Vancoppenolle (UCL-ASTR, Belgium) |
---|
| 12 | !! * Sylvain Bouillon (UCL-ASTR, Belgium) |
---|
| 13 | !! * Miguel Angel Morales Maqueda (NOC-L, UK) |
---|
| 14 | !! thanks to valuable earlier work by |
---|
| 15 | !! * Thierry Fichefet |
---|
| 16 | !! * Hugues Goosse |
---|
| 17 | !! thanks also to the following persons who contributed |
---|
| 18 | !! * Gurvan Madec, Claude Talandier, Christian Ethe (LOCEAN, France) |
---|
| 19 | !! * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany) |
---|
| 20 | !! * Bill Lipscomb (LANL), Cecilia Bitz (UWa) and Elisabeth Hunke (LANL), USA. |
---|
| 21 | !! |
---|
| 22 | !! SI3 has been made possible by a handful of persons who met as working group |
---|
| 23 | !! (from France, Belgium, UK and Italy) |
---|
| 24 | !! * Clement Rousset, Martin Vancoppenolle & Gurvan Madec (LOCEAN, France) |
---|
| 25 | !! * Matthieu Chevalier & David Salas (Meteo France, France) |
---|
| 26 | !! * Gilles Garric (Mercator Ocean, France) |
---|
| 27 | !! * Thierry Fichefet & Francois Massonnet (UCL, Belgium) |
---|
| 28 | !! * Ed Blockley & Jeff Ridley (Met Office, UK) |
---|
| 29 | !! * Danny Feltham & David Schroeder (CPOM, UK) |
---|
| 30 | !! * Yevgeny Aksenov (NOC, UK) |
---|
| 31 | !! * Paul Holland (BAS, UK) |
---|
| 32 | !! * Dorotea Iovino (CMCC, Italy) |
---|
| 33 | !!====================================================================== |
---|
| 34 | !! History : 4.0 ! 2018 (C. Rousset) Original code SI3 |
---|
| 35 | !!---------------------------------------------------------------------- |
---|
| 36 | #if defined key_si3 |
---|
| 37 | !!---------------------------------------------------------------------- |
---|
| 38 | !! 'key_si3' SI3 sea-ice model |
---|
| 39 | !!---------------------------------------------------------------------- |
---|
| 40 | !! ice_stp : sea-ice model time-stepping and update ocean SBC over ice-covered area |
---|
| 41 | !! ice_init : initialize sea-ice |
---|
| 42 | !!---------------------------------------------------------------------- |
---|
| 43 | USE oce ! ocean dynamics and tracers |
---|
| 44 | USE dom_oce ! ocean space and time domain |
---|
| 45 | USE c1d ! 1D vertical configuration |
---|
| 46 | USE ice ! sea-ice: variables |
---|
| 47 | USE ice1D ! sea-ice: thermodynamical 1D variables |
---|
| 48 | ! |
---|
| 49 | USE phycst ! Define parameters for the routines |
---|
| 50 | USE eosbn2 ! equation of state |
---|
| 51 | USE sbc_oce ! Surface boundary condition: ocean fields |
---|
| 52 | USE sbc_ice ! Surface boundary condition: ice fields |
---|
| 53 | ! |
---|
| 54 | USE icesbc ! sea-ice: Surface boundary conditions |
---|
| 55 | USE iceupdate ! sea-ice: sea surface boundary condition update |
---|
| 56 | USE icewri ! sea-ice: outputs |
---|
| 57 | ! |
---|
| 58 | ! |
---|
| 59 | USE in_out_manager ! I/O manager |
---|
| 60 | USE iom ! I/O manager library |
---|
| 61 | USE lib_mpp ! MPP library |
---|
| 62 | USE lib_fortran ! fortran utilities (glob_sum + no signed zero) |
---|
| 63 | USE timing ! Timing |
---|
| 64 | USE prtctl ! Print control |
---|
| 65 | |
---|
| 66 | IMPLICIT NONE |
---|
| 67 | PRIVATE |
---|
| 68 | |
---|
| 69 | PUBLIC ice_stp ! called by sbcmod.F90 |
---|
| 70 | PUBLIC ice_init ! called by sbcmod.F90 |
---|
| 71 | |
---|
| 72 | !! * Substitutions |
---|
| 73 | # include "do_loop_substitute.h90" |
---|
| 74 | !!---------------------------------------------------------------------- |
---|
| 75 | !! NEMO/ICE 4.0 , NEMO Consortium (2018) |
---|
[13686] | 76 | !! $Id: icestp.F90 13655 2020-10-21 14:15:13Z laurent $ |
---|
[13655] | 77 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
| 78 | !!---------------------------------------------------------------------- |
---|
| 79 | CONTAINS |
---|
| 80 | |
---|
| 81 | SUBROUTINE ice_stp( kt, Kbb, Kmm, ksbc ) |
---|
| 82 | !!--------------------------------------------------------------------- |
---|
| 83 | !! *** ROUTINE ice_stp *** |
---|
| 84 | !! |
---|
| 85 | !! ** Purpose : sea-ice model time-stepping and update ocean surface |
---|
| 86 | !! boundary condition over ice-covered area |
---|
| 87 | !! |
---|
| 88 | !! ** Method : ice model time stepping |
---|
| 89 | !! - call the ice dynamics routine |
---|
| 90 | !! - call the ice advection/diffusion routine |
---|
| 91 | !! - call the ice thermodynamics routine |
---|
| 92 | !! - call the routine that computes mass and |
---|
| 93 | !! heat fluxes at the ice/ocean interface |
---|
| 94 | !! - save the outputs |
---|
| 95 | !! - save the outputs for restart when necessary |
---|
| 96 | !! |
---|
| 97 | !! ** Action : - time evolution of the LIM sea-ice model |
---|
| 98 | !! - update all sbc variables below sea-ice: |
---|
| 99 | !! utau, vtau, taum, wndm, qns , qsr, emp , sfx |
---|
| 100 | !!--------------------------------------------------------------------- |
---|
| 101 | INTEGER, INTENT(in) :: kt ! ocean time step |
---|
| 102 | INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices |
---|
| 103 | INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk, or Pure Coupled) |
---|
| 104 | ! |
---|
| 105 | INTEGER :: jl ! dummy loop index |
---|
| 106 | !!---------------------------------------------------------------------- |
---|
| 107 | ! |
---|
| 108 | IF( ln_timing ) CALL timing_start('ice_stp') |
---|
| 109 | ! |
---|
| 110 | ! !-----------------------! |
---|
| 111 | IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! --- Ice time step --- ! |
---|
| 112 | ! !-----------------------! |
---|
| 113 | ! |
---|
| 114 | kt_ice = kt ! -- Ice model time step |
---|
| 115 | ! |
---|
| 116 | u_oce(:,:) = ssu_m(:,:) ! -- mean surface ocean current |
---|
| 117 | v_oce(:,:) = ssv_m(:,:) |
---|
| 118 | ! |
---|
| 119 | CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) ! -- freezing temperature [Kelvin] (set to rt0 over land) |
---|
| 120 | t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) |
---|
| 121 | ! |
---|
| 122 | ! |
---|
| 123 | !------------------------------------------------! |
---|
| 124 | ! --- Dynamical coupling with the atmosphere --- ! |
---|
| 125 | !------------------------------------------------! |
---|
| 126 | ! It provides the following fields used in sea ice model: |
---|
| 127 | ! utau_ice, vtau_ice = surface ice stress [N/m2] |
---|
| 128 | !------------------------------------------------! |
---|
[13686] | 129 | CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice ) |
---|
| 130 | ! |
---|
[13655] | 131 | !------------------------------------------------------! |
---|
| 132 | ! --- Thermodynamical coupling with the atmosphere --- ! |
---|
| 133 | !------------------------------------------------------! |
---|
| 134 | ! It provides the following fields used in sea ice model: |
---|
| 135 | ! emp_oce , emp_ice = E-P over ocean and sea ice [Kg/m2/s] |
---|
| 136 | ! sprecip = solid precipitation [Kg/m2/s] |
---|
| 137 | ! evap_ice = sublimation [Kg/m2/s] |
---|
| 138 | ! qsr_tot , qns_tot = solar & non solar heat flux (total) [W/m2] |
---|
| 139 | ! qsr_ice , qns_ice = solar & non solar heat flux over ice [W/m2] |
---|
| 140 | ! dqns_ice = non solar heat sensistivity [W/m2] |
---|
| 141 | ! qemp_oce, qemp_ice, = sensible heat (associated with evap & precip) [W/m2] |
---|
| 142 | ! qprec_ice, qevap_ice |
---|
| 143 | !------------------------------------------------------! |
---|
| 144 | CALL ice_sbc_flx( kt, ksbc ) |
---|
| 145 | !----------------------------! |
---|
| 146 | ! --- ice thermodynamics --- ! |
---|
| 147 | !----------------------------! |
---|
| 148 | |
---|
| 149 | fr_i (:,:) = at_i(:,:) !#LB... |
---|
| 150 | |
---|
| 151 | !!#LB: lines stolen from "ice_update_flx()@iceupdate.F90" |
---|
| 152 | !!============================================================= |
---|
| 153 | |
---|
| 154 | ! --- case we bypass ice thermodynamics --- ! |
---|
| 155 | !IF( .NOT. ln_icethd ) THEN ! we suppose ice is impermeable => ocean is isolated from atmosphere |
---|
| 156 | qt_atm_oi (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) |
---|
| 157 | qt_oce_ai (:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:) |
---|
| 158 | emp_ice (:,:) = 0._wp |
---|
| 159 | qemp_ice (:,:) = 0._wp |
---|
| 160 | qevap_ice (:,:,:) = 0._wp |
---|
| 161 | !ENDIF |
---|
| 162 | |
---|
| 163 | ! output all fluxes |
---|
| 164 | !------------------ |
---|
| 165 | |
---|
| 166 | ! --- mass fluxes [kg/m2/s] --- ! |
---|
| 167 | 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) |
---|
| 168 | 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) |
---|
| 169 | ! --- heat fluxes [W/m2] --- ! |
---|
| 170 | IF( iom_use('qsr_oce' ) ) CALL iom_put( 'qsr_oce' , qsr_oce * ( 1._wp - at_i_b ) ) ! solar flux at ocean surface |
---|
| 171 | 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 |
---|
| 172 | 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 |
---|
| 173 | IF( iom_use('qsr_ice' ) ) CALL iom_put( 'qsr_ice' , SUM( qsr_ice * a_i_b, dim=3 ) ) |
---|
| 174 | IF( iom_use('qt_oce' ) ) CALL iom_put( 'qt_oce' , ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce ) |
---|
| 175 | IF( iom_use('qt_ice' ) ) CALL iom_put( 'qt_ice' , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 ) + qemp_ice ) |
---|
| 176 | IF( iom_use('qemp_oce' ) ) CALL iom_put( 'qemp_oce' , qemp_oce ) ! Downward Heat Flux from E-P over ocean |
---|
| 177 | IF( iom_use('qemp_ice' ) ) CALL iom_put( 'qemp_ice' , qemp_ice ) ! Downward Heat Flux from E-P over ice |
---|
| 178 | |
---|
| 179 | |
---|
| 180 | !!============================================================= |
---|
| 181 | |
---|
| 182 | ! |
---|
| 183 | CALL ice_wri( kt ) ! -- Ice outputs |
---|
| 184 | ! |
---|
| 185 | !IF( lrst_ice ) CALL ice_rst_write( kt ) ! -- Ice restart file |
---|
| 186 | ! |
---|
| 187 | !IF( ln_icectl ) CALL ice_ctl( kt ) ! -- Control checks |
---|
| 188 | ! |
---|
| 189 | ENDIF ! End sea-ice time step only |
---|
| 190 | |
---|
| 191 | !-------------------------! |
---|
| 192 | ! --- Ocean time step --- ! |
---|
| 193 | !-------------------------! |
---|
| 194 | ! |
---|
| 195 | IF( ln_timing ) CALL timing_stop('ice_stp') |
---|
| 196 | ! |
---|
| 197 | END SUBROUTINE ice_stp |
---|
| 198 | |
---|
| 199 | |
---|
| 200 | SUBROUTINE ice_init( Kbb, Kmm, Kaa ) |
---|
| 201 | !!---------------------------------------------------------------------- |
---|
| 202 | !! *** ROUTINE ice_init *** |
---|
| 203 | !! |
---|
| 204 | !! ** purpose : Initialize sea-ice parameters |
---|
| 205 | !!---------------------------------------------------------------------- |
---|
| 206 | INTEGER, INTENT(in) :: Kbb, Kmm, Kaa |
---|
| 207 | ! |
---|
| 208 | INTEGER :: ierr |
---|
| 209 | !!---------------------------------------------------------------------- |
---|
| 210 | IF(lwp) WRITE(numout,*) |
---|
| 211 | IF(lwp) WRITE(numout,*) 'Sea Ice Model: SI3 (Sea Ice modelling Integrated Initiative)' |
---|
| 212 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' |
---|
| 213 | IF(lwp) WRITE(numout,*) |
---|
| 214 | IF(lwp) WRITE(numout,*) 'ice_init: Arrays allocation & Initialization of all routines & init state' |
---|
| 215 | IF(lwp) WRITE(numout,*) '~~~~~~~~' |
---|
| 216 | ! |
---|
| 217 | ! ! Load the reference and configuration namelist files and open namelist output file |
---|
| 218 | CALL load_nml( numnam_ice_ref, 'namelist_ice_ref', numout, lwm ) |
---|
| 219 | CALL load_nml( numnam_ice_cfg, 'namelist_ice_cfg', numout, lwm ) |
---|
| 220 | IF(lwm) CALL ctl_opn( numoni , 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 ) |
---|
| 221 | ! |
---|
| 222 | CALL par_init ! set some ice run parameters |
---|
| 223 | ! |
---|
| 224 | ! ! Allocate the ice arrays (sbc_ice already allocated in sbc_init) |
---|
| 225 | ierr = ice_alloc () ! ice variables |
---|
| 226 | ierr = ierr + sbc_ice_alloc () ! surface boundary conditions |
---|
| 227 | ierr = ierr + ice1D_alloc () ! thermodynamics |
---|
| 228 | ! |
---|
| 229 | CALL mpp_sum( 'icestp', ierr ) |
---|
| 230 | IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays') |
---|
| 231 | ! |
---|
| 232 | ! |
---|
| 233 | CALL ice_sbc_init ! set ice-ocean and ice-atm. coupling parameters |
---|
| 234 | ! |
---|
[13686] | 235 | fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction |
---|
| 236 | ! |
---|
[13655] | 237 | IF( ln_rstart ) CALL iom_close( numrir ) ! close input ice restart file |
---|
| 238 | ! |
---|
| 239 | END SUBROUTINE ice_init |
---|
| 240 | |
---|
| 241 | |
---|
| 242 | SUBROUTINE par_init |
---|
| 243 | !!------------------------------------------------------------------- |
---|
| 244 | !! *** ROUTINE par_init *** |
---|
| 245 | !! |
---|
| 246 | !! ** Purpose : Definition generic parameters for ice model |
---|
| 247 | !! |
---|
| 248 | !! ** Method : Read namelist and check the parameter |
---|
| 249 | !! values called at the first timestep (nit000) |
---|
| 250 | !! |
---|
| 251 | !! ** input : Namelist nampar |
---|
| 252 | !!------------------------------------------------------------------- |
---|
| 253 | INTEGER :: ios ! Local integer |
---|
| 254 | !! |
---|
| 255 | NAMELIST/nampar/ jpl, nlay_i, nlay_s, ln_virtual_itd, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s, & |
---|
| 256 | & cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir |
---|
| 257 | !!------------------------------------------------------------------- |
---|
| 258 | ! |
---|
| 259 | READ ( numnam_ice_ref, nampar, IOSTAT = ios, ERR = 901) |
---|
| 260 | 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampar in reference namelist' ) |
---|
| 261 | READ ( numnam_ice_cfg, nampar, IOSTAT = ios, ERR = 902 ) |
---|
| 262 | 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nampar in configuration namelist' ) |
---|
| 263 | IF(lwm) WRITE( numoni, nampar ) |
---|
| 264 | ! |
---|
| 265 | IF(lwp) THEN ! control print |
---|
| 266 | WRITE(numout,*) |
---|
| 267 | WRITE(numout,*) ' par_init: ice parameters shared among all the routines' |
---|
| 268 | WRITE(numout,*) ' ~~~~~~~~' |
---|
| 269 | WRITE(numout,*) ' Namelist nampar: ' |
---|
| 270 | WRITE(numout,*) ' number of ice categories jpl = ', jpl |
---|
| 271 | WRITE(numout,*) ' number of ice layers nlay_i = ', nlay_i |
---|
| 272 | WRITE(numout,*) ' number of snow layers nlay_s = ', nlay_s |
---|
| 273 | WRITE(numout,*) ' virtual ITD param for jpl=1 (T) or not (F) ln_virtual_itd = ', ln_virtual_itd |
---|
| 274 | WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_icedyn = ', ln_icedyn |
---|
| 275 | WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_icethd = ', ln_icethd |
---|
| 276 | WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n |
---|
| 277 | WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s |
---|
| 278 | ENDIF |
---|
| 279 | ! !--- change max ice concentration for roundoff errors |
---|
| 280 | rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 ) |
---|
| 281 | rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 ) |
---|
| 282 | ! !--- check consistency |
---|
| 283 | IF ( jpl > 1 .AND. ln_virtual_itd ) THEN |
---|
| 284 | ln_virtual_itd = .FALSE. |
---|
| 285 | IF(lwp) WRITE(numout,*) |
---|
| 286 | IF(lwp) WRITE(numout,*) ' ln_virtual_itd forced to false as jpl>1, no need with multiple categories to emulate them' |
---|
| 287 | ENDIF |
---|
| 288 | ! |
---|
| 289 | IF( ln_cpl .AND. nn_cats_cpl /= 1 .AND. nn_cats_cpl /= jpl ) THEN |
---|
| 290 | CALL ctl_stop( 'STOP', 'par_init: in coupled mode, nn_cats_cpl should be either 1 or jpl' ) |
---|
| 291 | ENDIF |
---|
| 292 | ! |
---|
| 293 | rDt_ice = REAL(nn_fsbc) * rn_Dt !--- sea-ice timestep and its inverse |
---|
| 294 | r1_Dt_ice = 1._wp / rDt_ice |
---|
| 295 | IF(lwp) WRITE(numout,*) |
---|
| 296 | IF(lwp) WRITE(numout,*) ' ice timestep rDt_ice = nn_fsbc*rn_Dt = ', rDt_ice |
---|
| 297 | ! |
---|
| 298 | r1_nlay_i = 1._wp / REAL( nlay_i, wp ) !--- inverse of nlay_i and nlay_s |
---|
| 299 | r1_nlay_s = 1._wp / REAL( nlay_s, wp ) |
---|
| 300 | ! |
---|
| 301 | END SUBROUTINE par_init |
---|
| 302 | |
---|
| 303 | #else |
---|
| 304 | !!---------------------------------------------------------------------- |
---|
| 305 | !! Default option Dummy module NO SI3 sea-ice model |
---|
| 306 | !!---------------------------------------------------------------------- |
---|
| 307 | CONTAINS |
---|
| 308 | SUBROUTINE ice_stp ( kt, ksbc ) ! Dummy routine |
---|
| 309 | INTEGER, INTENT(in) :: kt, ksbc |
---|
| 310 | WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt |
---|
| 311 | END SUBROUTINE ice_stp |
---|
| 312 | SUBROUTINE ice_init ! Dummy routine |
---|
| 313 | END SUBROUTINE ice_init |
---|
| 314 | #endif |
---|
| 315 | |
---|
| 316 | !!====================================================================== |
---|
| 317 | END MODULE icestp |
---|