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) |
---|
76 | !! $Id: icestp.F90 13655 2020-10-21 14:15:13Z laurent $ |
---|
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 | !------------------------------------------------! |
---|
129 | CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice ) |
---|
130 | ! |
---|
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 | ! |
---|
235 | fr_i (:,:) = at_i(:,:) ! initialisation of sea-ice fraction |
---|
236 | ! |
---|
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 |
---|