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 @ 13338

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

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

Download in other formats: