New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

iceistate.F90 in NEMO/branches/2020/ticket2487/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/ticket2487/src/ICE/iceistate.F90 @ 15410

Last change on this file since 15410 was 15410, checked in by smueller, 3 years ago
Synchronizing with ^{/NEMO/releases/r4.0/r4.0-HEAD@15405 (ticket #2487)}
Property svn:keywords set to `Id`
File size: 29.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: