Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

iceistate.F90

Last change on this file was 13442, checked in by dancopsey, 4 years ago
Put majority of code in to get melt pond freezing to feed back into ice thickness.
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	s_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	IF( ln_iceini ) THEN
174	! !---------------!
175	IF( nn_iceini_file == 1 )THEN ! Read a file !
176	! !---------------!
177	WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp
178	ELSEWHERE ; zswitch(:,:) = 0._wp
179	END WHERE
180	!
181	CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
182	!
183	! -- mandatory fields -- !
184	zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
185	zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
186	zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
187
188	! -- optional fields -- !
189	! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
190	!
191	! ice salinity
192	IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
193	& si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
194	zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
195	!
196	! ice temperature
197	IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) &
198	& si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
199	ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
200	!
201	! surface temperature => set to ice temperature if it exists
202	IF ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
203	si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
204	ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
205	si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
206	ENDIF
207	zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
208	!
209	! snow temperature => set to ice temperature if it exists
210	IF ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
211	si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
212	ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
213	si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
214	ENDIF
215	ztm_s_ini(:,:) = si(jp_tms)%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	zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
222	!
223	! pond depth
224	IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
225	& si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
226	zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
227	!
228	! change the switch for the following
229	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
230	ELSEWHERE ; zswitch(:,:) = 0._wp
231	END WHERE
232	! !---------------!
233	ELSE ! Read namelist !
234	! !---------------!
235	! no ice if (sst - Tfreez) >= thresold
236	WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
237	ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
238	END WHERE
239	!
240	! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
241	WHERE( ff_t(:,:) >= 0._wp )
242	zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
243	zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
244	zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
245	zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
246	ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
247	zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
248	ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
249	zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
250	zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
251	ELSEWHERE
252	zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
253	zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
254	zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
255	zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
256	ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
257	zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
258	ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
259	zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
260	zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
261	END WHERE
262	!
263	ENDIF
264
265	! make sure ponds = 0 if no ponds scheme
266	IF ( .NOT.ln_pnd ) THEN
267	zapnd_ini(:,:) = 0._wp
268	zhpnd_ini(:,:) = 0._wp
269	ENDIF
270
271	!-------------!
272	! fill fields !
273	!-------------!
274	! select ice covered grid points
275	npti = 0 ; nptidx(:) = 0
276	DO jj = 1, jpj
277	DO ji = 1, jpi
278	IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
279	npti = npti + 1
280	nptidx(npti) = (jj - 1) * jpi + ji
281	ENDIF
282	END DO
283	END DO
284
285	! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
286	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini )
287	CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini )
288	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini )
289	CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
290	CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
291	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini )
292	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini )
293	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini )
294	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini )
295
296	! allocate temporary arrays
297	ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
298	& zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
299
300	! distribute 1-cat into jpl-cat: (jpijpj) -> (jpijpj,jpl)
301	CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), &
302	& zhi_2d , zhs_2d , zai_2d , &
303	& 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), &
304	& zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d )
305
306	! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
307	DO jl = 1, jpl
308	zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
309	zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
310	END DO
311	CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i )
312	CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s )
313	CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i )
314	CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d )
315	CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d )
316	CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
317	CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i )
318	CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip )
319	CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip )
320
321	! deallocate temporary arrays
322	DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
323	& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
324
325	! calculate extensive and intensive variables
326	CALL ice_var_salprof ! for sz_i
327	DO jl = 1, jpl
328	DO jj = 1, jpj
329	DO ji = 1, jpi
330	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
331	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
332	sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
333	END DO
334	END DO
335	END DO
336	!
337	DO jl = 1, jpl
338	DO jk = 1, nlay_s
339	DO jj = 1, jpj
340	DO ji = 1, jpi
341	t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
342	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
343	& rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
344	END DO
345	END DO
346	END DO
347	END DO
348	!
349	DO jl = 1, jpl
350	DO jk = 1, nlay_i
351	DO jj = 1, jpj
352	DO ji = 1, jpi
353	t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
354	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
355	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
356	& rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
357	& rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
358	& - rcp * ( ztmelts - rt0 ) )
359	END DO
360	END DO
361	END DO
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	a_ip_eff(:,:,:) = a_ip_frac(:,:,:)
371	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
372
373	! specific temperatures for coupled runs
374	tn_ice(:,:,:) = t_su(:,:,:)
375	t1_ice(:,:,:) = t_i (:,:,1,:)
376	!
377	ENDIF ! ln_iceini
378	!
379	at_i(:,:) = SUM( a_i, dim=3 )
380	!
381	!----------------------------------------------
382	! 3) Snow-ice mass (case ice is fully embedded)
383	!----------------------------------------------
384	snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass
385	snwice_mass_b(:,:) = snwice_mass(:,:)
386	!
387	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
388	!
389	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
390	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
391	!
392	IF( .NOT.ln_linssh ) THEN
393	!
394	WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
395	ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE
396	!
397	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
398	e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
399	e3t_b(:,:,jk) = e3t_n(:,:,jk)
400	e3t_a(:,:,jk) = e3t_n(:,:,jk)
401	END DO
402	!
403	! Reconstruction of all vertical scale factors at now and before time-steps
404	! =========================================================================
405	! Horizontal scale factor interpolations
406	! --------------------------------------
407	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
408	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
409	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
410	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
411	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
412	! Vertical scale factor interpolations
413	! ------------------------------------
414	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
415	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
416	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
417	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
418	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
419	! t- and w- points depth
420	! ----------------------
421	!!gm not sure of that....
422	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
423	gdepw_n(:,:,1) = 0.0_wp
424	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
425	DO jk = 2, jpk
426	gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk )
427	gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
428	gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:)
429	END DO
430	ENDIF
431	ENDIF
432
433	!------------------------------------
434	! 4) store fields at before time-step
435	!------------------------------------
436	! it is only necessary for the 1st interpolation by Agrif
437	a_i_b (:,:,:) = a_i (:,:,:)
438	e_i_b (:,:,:,:) = e_i (:,:,:,:)
439	v_i_b (:,:,:) = v_i (:,:,:)
440	v_s_b (:,:,:) = v_s (:,:,:)
441	e_s_b (:,:,:,:) = e_s (:,:,:,:)
442	sv_i_b (:,:,:) = sv_i (:,:,:)
443	oa_i_b (:,:,:) = oa_i (:,:,:)
444	u_ice_b(:,:) = u_ice(:,:)
445	v_ice_b(:,:) = v_ice(:,:)
446	! total concentration is needed for Lupkes parameterizations
447	at_i_b (:,:) = at_i (:,:)
448
449	!!clem: output of initial state should be written here but it is impossible because
450	!! the ocean and ice are in the same file
451	!! CALL dia_wri_state( 'output.init' )
452	!
453	END SUBROUTINE ice_istate
454
455
456	SUBROUTINE ice_istate_init
457	!!-------------------------------------------------------------------
458	!! * ROUTINE ice_istate_init *
459	!!
460	!! ** Purpose : Definition of initial state of the ice
461	!!
462	!! ** Method : Read the namini namelist and check the parameter
463	!! values called at the first timestep (nit000)
464	!!
465	!! ** input : Namelist namini
466	!!
467	!!-----------------------------------------------------------------------------
468	INTEGER :: ios ! Local integer output status for namelist read
469	INTEGER :: ifpr, ierror
470	!
471	CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
472	TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
473	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
474	!
475	NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, &
476	& rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
477	& rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
478	& rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
479	& rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
480	& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
481	!!-----------------------------------------------------------------------------
482	!
483	REWIND( numnam_ice_ref ) ! Namelist namini in reference namelist : Ice initial state
484	READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
485	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' )
486	REWIND( numnam_ice_cfg ) ! Namelist namini in configuration namelist : Ice initial state
487	READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
488	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' )
489	IF(lwm) WRITE ( numoni, namini )
490	!
491	slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
492	slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi
493	slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms
494	slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd
495	!
496	IF(lwp) THEN ! control print
497	WRITE(numout,*)
498	WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
499	WRITE(numout,*) '~~~~~~~~~~~~~~~'
500	WRITE(numout,*) ' Namelist namini:'
501	WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini
502	WRITE(numout,*) ' ice initialization from a netcdf file nn_iceini_file = ', nn_iceini_file
503	WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
504	IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN
505	WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s
506	WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s
507	WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s
508	WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s
509	WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s
510	WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s
511	WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s
512	WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s
513	WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s
514	ENDIF
515	ENDIF
516	!
517	IF( nn_iceini_file == 1 ) THEN ! Ice initialization using input file
518	!
519	! set si structure
520	ALLOCATE( si(jpfldi), STAT=ierror )
521	IF( ierror > 0 ) THEN
522	CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN
523	ENDIF
524	!
525	DO ifpr = 1, jpfldi
526	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
527	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
528	END DO
529	!
530	! fill si with slf_i and control print
531	CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
532	!
533	ENDIF
534	!
535	IF( .NOT.ln_pnd ) THEN
536	rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
537	rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
538	CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
539	ENDIF
540	!
541	END SUBROUTINE ice_istate_init
542
543	#else
544	!!----------------------------------------------------------------------
545	!! Default option : Empty module NO SI3 sea-ice model
546	!!----------------------------------------------------------------------
547	#endif
548
549	!!======================================================================
550	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_retain_melt/src/ICE/iceistate.F90

Download in other formats: