Context Navigation

iceistate.F90 @ 12566

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/iceistate.F90 @ 12566

Download in other formats: