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/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/ICE/iceistate.F90 @ 13334

Last change on this file since 13334 was 13334, checked in by jchanut, 4 years ago
finish bypassing ocean/ice initialization with AGRIF, #2222, #2129
Property svn:keywords set to `Id`
File size: 28.5 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
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	USE agrif_oce ! initial state interpolation
36	USE agrif_ice_interp
37
38	IMPLICIT NONE
39	PRIVATE
40
41	PUBLIC ice_istate ! called by icestp.F90
42	PUBLIC ice_istate_init ! called by icestp.F90
43	!
44	! !! namelist (namini)
45	LOGICAL, PUBLIC :: ln_iceini !: Ice initialization or not
46	LOGICAL, PUBLIC :: ln_iceini_file !: Ice initialization from 2D netcdf 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
51	REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s
52	!
53	! ! if ln_iceini_file = T
54	INTEGER , PARAMETER :: jpfldi = 9 ! 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	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)
65
66	!! * Substitutions
67	# include "do_loop_substitute.h90"
68	!!----------------------------------------------------------------------
69	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
70	!! $Id$
71	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
72	!!----------------------------------------------------------------------
73	CONTAINS
74
75	SUBROUTINE ice_istate( kt, Kbb, Kmm, Kaa )
76	!!-------------------------------------------------------------------
77	!! * ROUTINE ice_istate *
78	!!
79	!! ** Purpose : defined the sea-ice initial state
80	!!
81	!! ** Method : This routine will put some ice where ocean
82	!! is at the freezing point, then fill in ice
83	!! state variables using prescribed initial
84	!! values in the namelist
85	!!
86	!! ** Steps : 1) Set initial surface and basal temperatures
87	!! 2) Recompute or read sea ice state variables
88	!! 3) Fill in the ice thickness distribution using gaussian
89	!! 4) Fill in space-dependent arrays for state variables
90	!! 5) snow-ice mass computation
91	!! 6) store before fields
92	!!
93	!! ** Notes : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
94	!! where there is no ice
95	!!--------------------------------------------------------------------
96	INTEGER, INTENT(in) :: kt ! time step
97	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
98	!
99	INTEGER :: ji, jj, jk, jl ! dummy loop indices
100	REAL(wp) :: ztmelts
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 !data from namelist or nc file
107	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !locak 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
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	a_ip_frac(:,:,:) = 0._wp
164	h_ip (:,:,:) = 0._wp
165	!
166	! ice velocities
167	u_ice (:,:) = 0._wp
168	v_ice (:,:) = 0._wp
169	!
170	!------------------------------------------------------------------------
171	! 2) overwrite some of the fields with namelist parameters or netcdf file
172	!------------------------------------------------------------------------
173
174
175	IF( ln_iceini ) THEN
176
177	#if defined key_agrif
178	IF ( ( Agrif_Root() ).OR.(.NOT.ln_init_chfrpar ) ) THEN
179	#endif
180	! !---------------!
181	IF( ln_iceini_file )THEN ! Read a file !
182	! !---------------!
183	WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp
184	ELSEWHERE ; zswitch(:,:) = 0._wp
185	END WHERE
186	!
187	CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
188	!
189	! -- mandatory fields -- !
190	zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
191	zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
192	zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)
193
194	! -- optional fields -- !
195	! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
196	!
197	! ice salinity
198	IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
199	& si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
200	!
201	! temperatures
202	IF ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. &
203	& TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN
204	si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
205	si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
206	si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
207	ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2
208	si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 )
209	ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2
210	si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 )
211	ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s
212	si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1)
213	ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i
214	si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
215	ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su
216	si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1)
217	ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i
218	si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
219	ENDIF
220	!
221	! pond concentration
222	IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
223	& 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.
224	& * si(jp_ati)%fnow(:,:,1)
225	!
226	! pond depth
227	IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
228	& si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_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	!
237	! change the switch for the following
238	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
239	ELSEWHERE ; zswitch(:,:) = 0._wp
240	END WHERE
241
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	ELSEWHERE
262	zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
263	zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
264	zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
265	zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
266	ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
267	zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
268	ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
269	zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
270	zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
271	END WHERE
272	!
273	ENDIF
274
275
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	ENDIF
282
283	!-------------!
284	! fill fields !
285	!-------------!
286	! select ice covered grid points
287	npti = 0 ; nptidx(:) = 0
288	DO_2D( 1, 1, 1, 1 )
289	IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
290	npti = npti + 1
291	nptidx(npti) = (jj - 1) * jpi + ji
292	ENDIF
293	END_2D
294
295	! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
296	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini )
297	CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini )
298	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini )
299	CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
300	CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
301	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini )
302	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini )
303	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini )
304	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini )
305
306	! allocate temporary arrays
307	ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
308	& zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
309
310	! distribute 1-cat into jpl-cat: (jpijpj) -> (jpijpj,jpl)
311	CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), &
312	& zhi_2d , zhs_2d , zai_2d , &
313	& 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), &
314	& zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d )
315
316	! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
317	DO jl = 1, jpl
318	zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
319	zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
320	END DO
321	CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i )
322	CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s )
323	CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i )
324	CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d )
325	CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d )
326	CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
327	CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i )
328	CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip )
329	CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip )
330
331	! deallocate temporary arrays
332	DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
333	& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
334
335	! calculate extensive and intensive variables
336	CALL ice_var_salprof ! for sz_i
337	DO jl = 1, jpl
338	DO_2D( 1, 1, 1, 1 )
339	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
340	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
341	sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
342	END_2D
343	END DO
344	!
345	DO jl = 1, jpl
346	DO_3D( 1, 1, 1, 1, 1, nlay_s )
347	t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
348	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
349	& rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
350	END_3D
351	END DO
352	!
353	DO jl = 1, jpl
354	DO_3D( 1, 1, 1, 1, 1, nlay_i )
355	t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
356	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
357	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
358	& rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
359	& rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
360	& - rcp * ( ztmelts - rt0 ) )
361	END_3D
362	END DO
363
364	! Melt ponds
365	WHERE( a_i > epsi10 )
366	a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
367	ELSEWHERE
368	a_ip_frac(:,:,:) = 0._wp
369	END WHERE
370	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
371
372	! specific temperatures for coupled runs
373	tn_ice(:,:,:) = t_su(:,:,:)
374	t1_ice(:,:,:) = t_i (:,:,1,:)
375	!
376	#if defined key_agrif
377	ELSE
378	CALL agrif_istate_ice
379	ENDIF
380	#endif
381	ENDIF ! ln_iceini
382	!
383	at_i(:,:) = SUM( a_i, dim=3 )
384	!
385	!----------------------------------------------
386	! 3) Snow-ice mass (case ice is fully embedded)
387	!----------------------------------------------
388	snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass
389	snwice_mass_b(:,:) = snwice_mass(:,:)
390	!
391	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
392	!
393	ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rho0
394	ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rho0
395	!
396	IF( .NOT.ln_linssh ) CALL dom_vvl_zgr( Kbb, Kmm, Kaa ) ! interpolation scale factor, depth and water column
397	! !!st
398	! IF( .NOT.ln_linssh ) THEN
399	! !
400	! WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:)
401	! ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE
402	! !
403	! DO jk = 1,jpkm1 ! adjust initial vertical scale factors
404	! e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:)
405	! e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
406	! e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm)
407	! END DO
408	! !
409	! ! Reconstruction of all vertical scale factors at now and before time-steps
410	! ! =========================================================================
411	! ! Horizontal scale factor interpolations
412	! ! --------------------------------------
413	! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
414	! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
415	! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
416	! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
417	! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
418	! ! Vertical scale factor interpolations
419	! ! ------------------------------------
420	! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' )
421	! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
422	! CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
423	! CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
424	! CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
425	! ! t- and w- points depth
426	! ! ----------------------
427	! !!gm not sure of that....
428	! gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
429	! gdepw(:,:,1,Kmm) = 0.0_wp
430	! gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
431	! DO jk = 2, jpk
432	! gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk ,Kmm)
433	! gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
434	! gde3w(:,:,jk) = gdept(:,:,jk ,Kmm) - ssh (:,:,Kmm)
435	! END DO
436	! ENDIF
437	ENDIF
438
439	!------------------------------------
440	! 4) store fields at before time-step
441	!------------------------------------
442	! it is only necessary for the 1st interpolation by Agrif
443	a_i_b (:,:,:) = a_i (:,:,:)
444	e_i_b (:,:,:,:) = e_i (:,:,:,:)
445	v_i_b (:,:,:) = v_i (:,:,:)
446	v_s_b (:,:,:) = v_s (:,:,:)
447	e_s_b (:,:,:,:) = e_s (:,:,:,:)
448	sv_i_b (:,:,:) = sv_i (:,:,:)
449	oa_i_b (:,:,:) = oa_i (:,:,:)
450	u_ice_b(:,:) = u_ice(:,:)
451	v_ice_b(:,:) = v_ice(:,:)
452	! total concentration is needed for Lupkes parameterizations
453	at_i_b (:,:) = at_i (:,:)
454
455	!!clem: output of initial state should be written here but it is impossible because
456	!! the ocean and ice are in the same file
457	!! CALL dia_wri_state( Kmm, 'output.init' )
458	!
459	END SUBROUTINE ice_istate
460
461
462	SUBROUTINE ice_istate_init
463	!!-------------------------------------------------------------------
464	!! * ROUTINE ice_istate_init *
465	!!
466	!! ** Purpose : Definition of initial state of the ice
467	!!
468	!! ** Method : Read the namini namelist and check the parameter
469	!! values called at the first timestep (nit000)
470	!!
471	!! ** input : Namelist namini
472	!!
473	!!-----------------------------------------------------------------------------
474	INTEGER :: ios, ifpr, ierror ! Local integers
475
476	!
477	CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
478	TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
479	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
480	!
481	NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
482	& rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
483	& rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
484	& rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
485	& rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
486	& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
487	!!-----------------------------------------------------------------------------
488	!
489	READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
490	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' )
491	READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
492	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' )
493	IF(lwm) WRITE ( numoni, namini )
494	!
495	slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
496	slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi
497	slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms
498	slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd
499	!
500	IF(lwp) THEN ! control print
501	WRITE(numout,*)
502	WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
503	WRITE(numout,*) '~~~~~~~~~~~~~~~'
504	WRITE(numout,*) ' Namelist namini:'
505	WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini
506	WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file
507	WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
508	IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
509	WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s
510	WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s
511	WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s
512	WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s
513	WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s
514	WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s
515	WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s
516	WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s
517	WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s
518	ENDIF
519	ENDIF
520	!
521	IF( ln_iceini_file ) THEN ! Ice initialization using input file
522	!
523	! set si structure
524	ALLOCATE( si(jpfldi), STAT=ierror )
525	IF( ierror > 0 ) THEN
526	CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN
527	ENDIF
528	!
529	DO ifpr = 1, jpfldi
530	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
531	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
532	END DO
533	!
534	! fill si with slf_i and control print
535	CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
536	!
537	ENDIF
538	!
539	IF( .NOT.ln_pnd ) THEN
540	rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
541	rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
542	CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
543	ENDIF
544	!
545	END SUBROUTINE ice_istate_init
546
547	#else
548	!!----------------------------------------------------------------------
549	!! Default option : Empty module NO SI3 sea-ice model
550	!!----------------------------------------------------------------------
551	#endif
552
553	!!======================================================================
554	END MODULE iceistate

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

Download in other formats: