Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
MODULE iceistate
!!======================================================================
!! *** MODULE iceistate ***
!! sea-ice : Initialization of ice variables
!!======================================================================
!! History : 2.0 ! 2004-01 (C. Ethe, G. Madec) Original code
!! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats
!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!----------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_istate : initialization of diagnostics ice variables
!! ice_istate_init : initialization of ice state and namelist read
!!----------------------------------------------------------------------
USE phycst ! physical constant
USE oce ! dynamics and tracers variables
USE dom_oce ! ocean domain
USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd
USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
USE eosbn2 ! equation of state
# if defined key_qco
USE domqco ! Quasi-Eulerian coord.
# elif defined key_linssh
! ! Fix in time coord.
# else
USE domvvl ! Variable volume
# endif
USE ice ! sea-ice: variables
USE ice1D ! sea-ice: thermodynamics variables
USE icetab ! sea-ice: 1D <==> 2D transformation
USE icevar ! sea-ice: operations
!
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 fldread ! read input fields
# if defined key_agrif
USE agrif_oce
USE agrif_ice_interp
# endif
IMPLICIT NONE
PRIVATE
PUBLIC ice_istate ! called by icestp.F90
PUBLIC ice_istate_init ! called by icestp.F90
!
! !! ** namelist (namini) **
LOGICAL, PUBLIC :: ln_iceini !: Ice initialization or not
INTEGER, PUBLIC :: nn_iceini_file !: Ice initialization:
! 0 = Initialise sea ice based on SSTs
! 1 = Initialise sea ice from single category netcdf file
! 2 = Initialise sea ice from multi category restart file
REAL(wp) :: rn_thres_sst
REAL(wp) :: rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
REAL(wp) :: rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
REAL(wp) :: rn_apd_ini_n, rn_hpd_ini_n, rn_hld_ini_n
REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s, rn_hld_ini_s
!
! ! if nn_iceini_file = 1
INTEGER , PARAMETER :: jpfldi = 10 ! maximum number of files to read
INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m)
INTEGER , PARAMETER :: jp_hts = 2 ! index of snw thickness (m)
INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (-)
INTEGER , PARAMETER :: jp_smi = 4 ! index of ice salinity (g/kg)
INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temperature (K)
INTEGER , PARAMETER :: jp_tsu = 6 ! index of ice surface temp (K)
INTEGER , PARAMETER :: jp_tms = 7 ! index of snw temperature (K)
INTEGER , PARAMETER :: jp_apd = 8 ! index of pnd fraction (-)
INTEGER , PARAMETER :: jp_hpd = 9 ! index of pnd depth (m)
INTEGER , PARAMETER :: jp_hld = 10 ! index of pnd lid depth (m)
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)
!
#if defined key_agrif
REAL(wp), PUBLIC :: rsshadj !: initial mean ssh adjustment due to initial ice+snow mass
#endif
!
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: iceistate.F90 15530 2021-11-23 15:09:32Z clem $
!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_istate( kt, Kbb, Kmm, Kaa )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_istate ***
!!
!! ** Purpose : defined the sea-ice initial state
!!
!! ** Method : This routine will put some ice where ocean
!! is at the freezing point, then fill in ice
!! state variables using prescribed initial
!! values in the namelist
!!
!! ** Steps : 1) Set initial surface and basal temperatures
!! 2) Recompute or read sea ice state variables
!! 3) Fill in space-dependent arrays for state variables
!! 4) snow-ice mass computation
!!
!! ** Notes : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
!! where there is no ice
!!--------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! time step
INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
!
INTEGER :: ji, jj, jk, jl ! dummy loop indices
REAL(wp) :: ztmelts, zsshadj, area
INTEGER , DIMENSION(4) :: itest
REAL(wp), DIMENSION(jpi,jpj) :: z2d
REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator
REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file
REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini, zhlid_ini !data from namelist or nc file
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays
!!
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d
!--------------------------------------------------------------------
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
!---------------------------
! 1) 1st init. of the fields
!---------------------------
!
! basal temperature (considered at freezing point) [Kelvin]
CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
!
! surface temperature and conductivity
DO jl = 1, jpl
t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface
cnd_ice(:,:,jl) = 0._wp ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
END DO
!
! ice and snw temperatures
DO jl = 1, jpl
DO jk = 1, nlay_i
t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
END DO
DO jk = 1, nlay_s
t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
END DO
END DO
!
! specific temperatures for coupled runs
tn_ice (:,:,:) = t_i (:,:,1,:)
t1_ice (:,:,:) = t_i (:,:,1,:)
! heat contents
e_i (:,:,:,:) = 0._wp
e_s (:,:,:,:) = 0._wp
! general fields
a_i (:,:,:) = 0._wp
v_i (:,:,:) = 0._wp
v_s (:,:,:) = 0._wp
sv_i(:,:,:) = 0._wp
oa_i(:,:,:) = 0._wp
!
h_i (:,:,:) = 0._wp
h_s (:,:,:) = 0._wp
s_i (:,:,:) = 0._wp
o_i (:,:,:) = 0._wp
!
! melt ponds
a_ip (:,:,:) = 0._wp
v_ip (:,:,:) = 0._wp
v_il (:,:,:) = 0._wp
a_ip_eff (:,:,:) = 0._wp
h_ip (:,:,:) = 0._wp
h_il (:,:,:) = 0._wp
!
! ice velocities
u_ice (:,:) = 0._wp
v_ice (:,:) = 0._wp
!
!------------------------------------------------------------------------
! 2) overwrite some of the fields with namelist parameters or netcdf file
!------------------------------------------------------------------------
IF( ln_iceini ) THEN
!
#if defined key_agrif
IF ( ( Agrif_Root() ).OR.(.NOT.ln_init_chfrpar ) ) THEN
#endif
! !---------------!
IF( nn_iceini_file == 1 )THEN ! Read a file !
! !---------------!
WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp
ELSEWHERE ; zswitch(:,:) = 0._wp
END WHERE
!
CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
!
! -- mandatory fields -- !
zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)
! -- optional fields -- !
! if fields do not exist then set them to the values present in the namelist (except for temperatures)
!
! ice salinity
IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
& si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
!
! temperatures
IF ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. &
& TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN
si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
ENDIF
IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
& si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
& si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_su, set T_su = T_s
& si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_su, set T_su = T_i
& si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_s, set T_s = T_su
& si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_s, set T_s = T_i
& si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
!
! pond concentration
IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
& si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
& * si(jp_ati)%fnow(:,:,1)
!
! pond depth
IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
& si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
!
! pond lid depth
IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) &
& si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
!
zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1)
ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1)
zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1)
ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1)
zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1)
zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1)
zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * tmask(:,:,1)
!
! change the switch for the following
WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
ELSEWHERE ; zswitch(:,:) = 0._wp
END WHERE
! !---------------!
ELSE ! Read namelist !
! !---------------!
! no ice if (sst - Tfreez) >= thresold
WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
END WHERE
!
! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
WHERE( ff_t(:,:) >= 0._wp )
zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:)
ELSEWHERE
zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:)
END WHERE
!
ENDIF
! make sure ponds = 0 if no ponds scheme
IF ( .NOT.ln_pnd ) THEN
zapnd_ini(:,:) = 0._wp
zhpnd_ini(:,:) = 0._wp
zhlid_ini(:,:) = 0._wp
ENDIF
IF ( .NOT.ln_pnd_lids ) THEN
zhlid_ini(:,:) = 0._wp
ENDIF
!----------------!
! 3) fill fields !
!----------------!
! select ice covered grid points
npti = 0 ; nptidx(:) = 0
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
npti = npti + 1
nptidx(npti) = (jj - 1) * jpi + ji
ENDIF
END_2D
! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini )
CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti) , zhlid_ini )
! allocate temporary arrays
ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), &
& zti_2d (npti,jpl), zts_2d (npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), &
& zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_2d(npti,jpl) )
! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), &
& zhi_2d , zhs_2d , zai_2d , &
& t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), &
& s_i_1d(1:npti) , a_ip_1d(1:npti) , h_ip_1d(1:npti), h_il_1d(1:npti), &
& zti_2d , zts_2d , ztsu_2d , &
& zsi_2d , zaip_2d , zhip_2d , zhil_2d )
! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
DO jl = 1, jpl
zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
END DO
CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i )
CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s )
CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i )
CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d )
CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d )
CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i )
CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip )
CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip )
CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d , h_il )
! deallocate temporary arrays
DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d )
! calculate extensive and intensive variables
CALL ice_var_salprof ! for sz_i
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
END_2D
END DO
!
DO jl = 1, jpl
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_s )
t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
& rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
END_3D
END DO
!
DO jl = 1, jpl
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
& rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
& rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
& - rcp * ( ztmelts - rt0 ) )
END_3D
END DO
!
#if defined key_agrif
ELSE
CALL agrif_istate_ice
ENDIF
#endif
! Melt ponds
WHERE( a_i > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
ELSEWHERE ; a_ip_eff(:,:,:) = 0._wp
END WHERE
v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
! specific temperatures for coupled runs
tn_ice(:,:,:) = t_su(:,:,:)
t1_ice(:,:,:) = t_i (:,:,1,:)
!
! ice concentration should not exceed amax
at_i(:,:) = SUM( a_i, dim=3 )
DO jl = 1, jpl
WHERE( at_i(:,:) > rn_amax_2d(:,:) ) a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:)
END DO
at_i(:,:) = SUM( a_i, dim=3 )
!
ENDIF ! ln_iceini
!
!----------------------------------------------------------
! 4) Adjust ssh and vertical scale factors to snow-ice mass
!----------------------------------------------------------
snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3 ) ! snow+ice mass
snwice_mass_b(:,:) = snwice_mass(:,:)
!
IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
! ! ----------------
ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rho0
ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0
!
ELSE ! levitating sea-ice: deplete the initial ssh over the whole domain
! ! ------------------
area = glob_sum( 'iceistate', e1e2t(:,:) * ssmask(:,:) )
zsshadj = glob_sum( 'iceistate', snwice_mass(:,:) * r1_rho0 * e1e2t(:,:) ) / area
#if defined key_agrif
! Override ssh adjustment in nested domains by the root-domain ssh adjustment;
! store the adjustment value in a global module variable to make it retrievable in nested domains
IF( .NOT.Agrif_Root() ) THEN
ELSE
rsshadj = zsshadj
ENDIF
#endif
IF(lwp) WRITE(numout,'(A23,F10.6,A20)') ' sea level adjusted by ', -zsshadj, ' m to compensate for'
IF(lwp) WRITE(numout,*) ' the initial snow+ice mass'
!
WHERE( ssmask(:,:) == 1._wp )
ssh(:,:,Kmm) = ssh(:,:,Kmm) - zsshadj
ssh(:,:,Kbb) = ssh(:,:,Kbb) - zsshadj
ENDWHERE
!
ENDIF
!
IF( .NOT.ln_linssh ) THEN
#if defined key_qco
CALL dom_qco_zgr( Kbb, Kmm ) ! upadte of r3=ssh/h0 ratios
#elif defined key_linssh
! ! Fix in time : key_linssh case, set through domzgr_substitute.h90
#else
DO jk = 1, jpk
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls)
e3t(ji,jj,jk,Kmm) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kmm) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
e3t(ji,jj,jk,Kbb) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kbb) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) )
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
END_2D
END DO
!
CALL dom_vvl_zgr( Kbb, Kmm, Kaa ) ! interpolation of scale factor, depth and water column
#endif
ENDIF
!
!!clem: output of initial state should be written here but it is impossible because
!! the ocean and ice are in the same file
!! CALL dia_wri_state( 'output.init' )
!
END SUBROUTINE ice_istate
SUBROUTINE ice_istate_init
!!-------------------------------------------------------------------
!! *** ROUTINE ice_istate_init ***
!!
!! ** Purpose : Definition of initial state of the ice
!!
!! ** Method : Read the namini namelist and check the parameter
!! values called at the first timestep (nit000)
!!
!! ** input : Namelist namini
!!
!!-----------------------------------------------------------------------------
INTEGER :: ios ! Local integer output status for namelist read
INTEGER :: ifpr, ierror
!
CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd, sn_hld
TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
!
NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, &
& rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
& rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
& rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
& rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, rn_hld_ini_n, rn_hld_ini_s, &
& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, sn_hld, cn_dir
!!-----------------------------------------------------------------------------
!
READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' )
READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' )
IF(lwm) WRITE ( numoni, namini )
!
slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi
slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms
slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd ; slf_i(jp_hld) = sn_hld
!
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
WRITE(numout,*) '~~~~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namini:'
WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini
WRITE(numout,*) ' ice initialization from a netcdf file nn_iceini_file = ', nn_iceini_file
WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN
WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s
WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s
WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s
WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s
WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s
WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s
WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s
WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s
WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s
WRITE(numout,*) ' initial pnd lid depth in the north-south rn_hld_ini = ', rn_hld_ini_n,rn_hld_ini_s
ENDIF
ENDIF
!
IF( nn_iceini_file == 1 ) THEN ! Ice initialization using input file
!
! set si structure
ALLOCATE( si(jpfldi), STAT=ierror )
IF( ierror > 0 ) THEN
CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN
ENDIF
!
DO ifpr = 1, jpfldi
ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
END DO
!
! fill si with slf_i and control print
CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
!
ENDIF
!
IF( .NOT.ln_pnd ) THEN
rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 & rn_hld_ini = 0 when no ponds' )
ENDIF
!
IF( .NOT.ln_pnd_lids ) THEN
rn_hld_ini_n = 0. ; rn_hld_ini_s = 0.
ENDIF
!
END SUBROUTINE ice_istate_init
#else
!!----------------------------------------------------------------------
!! Default option : Empty module NO SI3 sea-ice model
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE iceistate