Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

iceistate.F90 @ 13792

Last change on this file since 13792 was 12812, checked in by dancopsey, 4 years ago

Merge in Clem's branch:

http://forge.ipsl.jussieu.fr/nemo/browser/NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl

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

Download in other formats: