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
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
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
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
# 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) :: zswitch ! ice indicator
REAL(wp), DIMENSION(A2D(0)) :: zmsk ! ice indicator
REAL(wp), DIMENSION(A2D(0)) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file
REAL(wp), DIMENSION(A2D(0)) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
REAL(wp), DIMENSION(A2D(0)) :: 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)
!
DO jl = 1, jpl
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
! == reduced arrays == !
DO_2D( 0, 0, 0, 0 )
!
cnd_ice(ji,jj,jl) = 0._wp ! conductivity at the ice top
!
tn_ice(ji,jj,jl) = t_i (ji,jj,1,jl) ! temp for coupled runs
t1_ice(ji,jj,jl) = t_i (ji,jj,1,jl) ! temp for coupled runs
!
a_ip_eff(ji,jj,jl) = 0._wp ! melt pond effective fraction
END_2D
!
! == full arrays == !
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
! heat contents
DO jk = 1, nlay_i
e_i(ji,jj,jk,jl) = 0._wp
t_i(ji,jj,jk,jl) = rt0 * tmask(ji,jj,1) ! ice temp
ENDDO
DO jk = 1, nlay_s
e_s(ji,jj,jk,jl) = 0._wp
t_s(ji,jj,jk,jl) = rt0 * tmask(ji,jj,1) ! snw temp
ENDDO
!
! general fields
a_i (ji,jj,jl) = 0._wp
v_i (ji,jj,jl) = 0._wp
v_s (ji,jj,jl) = 0._wp
sv_i(ji,jj,jl) = 0._wp
oa_i(ji,jj,jl) = 0._wp
h_i (ji,jj,jl) = 0._wp
h_s (ji,jj,jl) = 0._wp
s_i (ji,jj,jl) = 0._wp
o_i (ji,jj,jl) = 0._wp
t_su(ji,jj,jl) = rt0 * tmask(ji,jj,1)
!
! melt ponds
a_ip(ji,jj,jl) = 0._wp
v_ip(ji,jj,jl) = 0._wp
v_il(ji,jj,jl) = 0._wp
h_ip(ji,jj,jl) = 0._wp
h_il(ji,jj,jl) = 0._wp
!
END_2D
!
ENDDO
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
u_ice(ji,jj) = 0._wp
v_ice(ji,jj) = 0._wp
END_2D
!
!------------------------------------------------------------------------
! 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(A2D(0)) >= 0._wp ) ; zmsk(:,:) = 1._wp
ELSEWHERE ; zmsk(:,:) = 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) * smask0(:,:)
zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * smask0(:,:)
zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * smask0(:,:)
! -- 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 * zmsk + rn_smi_ini_s * (1._wp - zmsk) ) * smask0(:,:)
!
! 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 * zmsk + rn_tmi_ini_s * (1._wp - zmsk) ) * smask0(:,:)
si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zmsk + rn_tsu_ini_s * (1._wp - zmsk) ) * smask0(:,:)
si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zmsk + rn_tms_ini_s * (1._wp - zmsk) ) * smask0(:,:)
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 * zmsk + rn_apd_ini_s * (1._wp - zmsk) ) * smask0(:,:) & ! 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 * zmsk + rn_hpd_ini_s * (1._wp - zmsk) ) * smask0(:,:)
!
! pond lid depth
IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) &
& si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zmsk + rn_hld_ini_s * (1._wp - zmsk) ) * smask0(:,:)
zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * smask0(:,:)
ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * smask0(:,:)
zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * smask0(:,:)
ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * smask0(:,:)
zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * smask0(:,:)
zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * smask0(:,:)
zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * smask0(:,:)
!
! !---------------!
ELSE ! Read namelist !
! !---------------!
! no ice if (sst - Tfreez) >= thresold
WHERE( ( sst_m(A2D(0)) - (t_bo(A2D(0)) - rt0) ) * smask0(:,:) >= rn_thres_sst ) ; zmsk(:,:) = 0._wp
ELSEWHERE ; zmsk(:,:) = smask0(:,:)
END WHERE
!
! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
WHERE( ff_t(A2D(0)) >= 0._wp )
zht_i_ini(:,:) = rn_hti_ini_n * zmsk(:,:)
zht_s_ini(:,:) = rn_hts_ini_n * zmsk(:,:)
zat_i_ini(:,:) = rn_ati_ini_n * zmsk(:,:)
zsm_i_ini(:,:) = rn_smi_ini_n * zmsk(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_n * zmsk(:,:)
zt_su_ini(:,:) = rn_tsu_ini_n * zmsk(:,:)
ztm_s_ini(:,:) = rn_tms_ini_n * zmsk(:,:)
zapnd_ini(:,:) = rn_apd_ini_n * zmsk(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
zhpnd_ini(:,:) = rn_hpd_ini_n * zmsk(:,:)
zhlid_ini(:,:) = rn_hld_ini_n * zmsk(:,:)
zht_i_ini(:,:) = rn_hti_ini_s * zmsk(:,:)
zht_s_ini(:,:) = rn_hts_ini_s * zmsk(:,:)
zat_i_ini(:,:) = rn_ati_ini_s * zmsk(:,:)
zsm_i_ini(:,:) = rn_smi_ini_s * zmsk(:,:)
ztm_i_ini(:,:) = rn_tmi_ini_s * zmsk(:,:)
zt_su_ini(:,:) = rn_tsu_ini_s * zmsk(:,:)
ztm_s_ini(:,:) = rn_tms_ini_s * zmsk(:,:)
zapnd_ini(:,:) = rn_apd_ini_s * zmsk(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
zhpnd_ini(:,:) = rn_hpd_ini_s * zmsk(:,:)
zhlid_ini(:,:) = rn_hld_ini_s * zmsk(:,:)
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
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
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 )
! this call is needed because of the calculations above that are done only in the interior
CALL lbc_lnk( 'iceistate', a_i , 'T', 1._wp, h_i , 'T', 1._wp, h_s , 'T', 1._wp, &
& zti_3d, 'T', 1._wp, zts_3d, 'T', 1._wp, t_su, 'T', 1._wp, &
& s_i , 'T', 1._wp, a_ip , 'T', 1._wp, h_ip, 'T', 1._wp, h_il, 'T', 1._wp )
! switch for the following
WHERE( SUM(a_i(:,:,:),dim=3) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
ELSEWHERE ; zswitch(:,:) = 0._wp
END WHERE
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
! 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(A2D(0),:) > epsi10 ) ; a_ip_eff(:,:,:) = a_ip(A2D(0),:) / a_i(A2D(0),:)
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
DO_2D( 0, 0, 0, 0 )
tn_ice(ji,jj,:) = t_su(ji,jj,:)
t1_ice(ji,jj,:) = t_i (ji,jj,1,:)
END_2D
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
!
! 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
zsshadj = Agrif_Parent(rsshadj)
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)
! IF( snwice_mass(ji,jj) /= 0._wp ) THEN
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) )
! ENDIF
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
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(A2D(0),1) )
IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(A2D(0),1,2) )
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
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