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.

icestp.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90 @ 8378

Last change on this file since 8378 was 8378, checked in by clem, 7 years ago
some cleaning
File size: 35.0 KB

Line
1	MODULE icestp
2	!!======================================================================
3	!! * MODULE icestp *
4	!! Surface module : update the ocean surface boundary condition over ice
5	!! & covered area using LIM sea-ice model
6	!! Sea-Ice model : LIM-3 Sea ice model time-stepping
7	!!=====================================================================
8	!! History : 2.0 ! 2006-12 (M. Vancoppenolle) Original code
9	!! 3.0 ! 2008-02 (C. Talandier) Surface module from icestp.F90
10	!! - ! 2008-04 (G. Madec) sltyle and lim_ctl routine
11	!! 3.3 ! 2010-11 (G. Madec) ice-ocean stress always computed at each ocean time-step
12	!! 3.4 ! 2011-01 (A Porter) dynamical allocation
13	!! - ! 2012-10 (C. Rousset) add lim_diahsb
14	!! 3.6 ! 2014-07 (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface
15	!! 4.0 ! 2016-06 (L. Brodeau) new unified bulk routine (based on AeroBulk)
16	!!----------------------------------------------------------------------
17	#if defined key_lim3
18	!!----------------------------------------------------------------------
19	!! 'key_lim3' : LIM 3.0 sea-ice model
20	!!----------------------------------------------------------------------
21	!! ice_stp : sea-ice model time-stepping and update ocean sbc over ice-covered area
22	!!----------------------------------------------------------------------
23	USE oce ! ocean dynamics and tracers
24	USE dom_oce ! ocean space and time domain
25	USE ice ! LIM-3: ice variables
26	USE thd_ice ! LIM-3: thermodynamical variables
27	!
28	USE sbc_oce ! Surface boundary condition: ocean fields
29	USE sbc_ice ! Surface boundary condition: ice fields
30	USE usrdef_sbc ! user defined: surface boundary condition
31	USE sbcblk ! Surface boundary condition: bulk
32	USE sbccpl ! Surface boundary condition: coupled interface
33	USE albedoice ! ice albedo
34	!
35	USE phycst ! Define parameters for the routines
36	USE eosbn2 ! equation of state
37	USE limdyn ! Ice dynamics
38	USE limtrp ! Ice transport
39	USE limthd ! Ice thermodynamics
40	USE limitd_me ! Mechanics on ice thickness distribution
41	USE limsbc ! sea surface boundary condition
42	USE limdiahsb ! Ice budget diagnostics
43	USE limwri ! Ice outputs
44	USE limrst ! Ice restarts
45	USE limupdate1 ! update of global variables
46	USE limupdate2 ! update of global variables
47	USE limvar ! Ice variables switch
48	USE limctl !
49	! MV MP 2016
50	USE limmp
51	! END MV MP 2016
52	USE limistate ! LIM initial state
53	USE limthd_sal ! LIM ice thermodynamics: salinity
54	!
55	USE c1d ! 1D vertical configuration
56	USE in_out_manager ! I/O manager
57	USE iom ! I/O manager library
58	USE prtctl ! Print control
59	USE lib_fortran !
60	USE lbclnk ! lateral boundary condition - MPP link
61	USE lib_mpp ! MPP library
62	USE timing ! Timing
63
64	USE bdy_oce , ONLY: ln_bdy
65	USE bdyice_lim ! unstructured open boundary data (bdy_ice_lim routine)
66	# if defined key_agrif
67	USE agrif_ice
68	USE agrif_lim3_update
69	USE agrif_lim3_interp
70	# endif
71
72	IMPLICIT NONE
73	PRIVATE
74
75	PUBLIC ice_stp ! routine called by sbcmod.F90
76	PUBLIC ice_init ! routine called by sbcmod.F90
77
78	!! * Substitutions
79	# include "vectopt_loop_substitute.h90"
80	!!----------------------------------------------------------------------
81	!! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
82	!! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $
83	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
84	!!----------------------------------------------------------------------
85	CONTAINS
86
87	SUBROUTINE ice_stp( kt, ksbc )
88	!!---------------------------------------------------------------------
89	!! * ROUTINE ice_stp *
90	!!
91	!! ** Purpose : update the ocean surface boundary condition via the
92	!! Louvain la Neuve Sea Ice Model time stepping
93	!!
94	!! ** Method : ice model time stepping
95	!! - call the ice dynamics routine
96	!! - call the ice advection/diffusion routine
97	!! - call the ice thermodynamics routine
98	!! - call the routine that computes mass and
99	!! heat fluxes at the ice/ocean interface
100	!! - save the outputs
101	!! - save the outputs for restart when necessary
102	!!
103	!! ** Action : - time evolution of the LIM sea-ice model
104	!! - update all sbc variables below sea-ice:
105	!! utau, vtau, taum, wndm, qns , qsr, emp , sfx
106	!!---------------------------------------------------------------------
107	INTEGER, INTENT(in) :: kt ! ocean time step
108	INTEGER, INTENT(in) :: ksbc ! type of sbc flux ( 1 = user defined formulation,
109	! 3 = bulk formulation,
110	! 4 = Pure Coupled formulation)
111	!!
112	INTEGER :: jl ! dummy loop index
113	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb_os, zalb_cs ! ice albedo under overcast/clear sky
114	REAL(wp), DIMENSION(jpi,jpj) :: zutau_ice, zvtau_ice
115	!!----------------------------------------------------------------------
116
117	IF( nn_timing == 1 ) CALL timing_start('ice_stp')
118
119	!-----------------------!
120	! --- Ice time step --- !
121	!-----------------------!
122	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
123
124	! Ice model time step
125	kt_ice = kt
126
127	# if defined key_agrif
128	IF( .NOT. Agrif_Root() ) lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
129	# endif
130
131	! mean surface ocean current at ice velocity point (C-grid dynamics : U- & V-points as the ocean)
132	u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
133	v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
134
135	! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
136	CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
137	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
138
139	! Mask sea ice surface temperature (set to rt0 over land)
140	DO jl = 1, jpl
141	t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
142	END DO
143	!
144	!------------------------------------------------!
145	! --- Dynamical coupling with the atmosphere --- !
146	!------------------------------------------------!
147	! It provides the following fields:
148	! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
149	!-----------------------------------------------------------------
150	CALL ice_bef ! Store previous ice values
151	SELECT CASE( ksbc )
152	CASE( jp_usr ) ; CALL usrdef_sbc_ice_tau( kt ) ! user defined formulation
153	CASE( jp_blk ) ; CALL blk_ice_tau ! Bulk formulation
154	CASE( jp_purecpl ) ; CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) ! Coupled formulation
155	END SELECT
156
157	IF( ln_mixcpl) THEN ! Case of a mixed Bulk/Coupled formulation
158	CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
159	utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
160	vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
161	ENDIF
162
163	!-------------------------------------------------------!
164	! --- ice dynamics and transport (except in 1D case) ---!
165	!-------------------------------------------------------!
166	CALL ice_diag0 ! set diag of mass, heat and salt fluxes to 0
167	CALL lim_rst_opn( kt ) ! Open Ice restart file
168	!
169	! --- zap this if no ice dynamics --- !
170	IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN
171	!
172	IF( nn_limdyn /= 0 ) THEN ! -- Ice dynamics
173	CALL lim_dyn( kt ) ! rheology
174	ELSE
175	u_ice(:,:) = rn_uice * umask(:,:,1) ! or prescribed velocity
176	v_ice(:,:) = rn_vice * vmask(:,:,1)
177	!!CALL RANDOM_NUMBER(u_ice(:,:))
178	!!CALL RANDOM_NUMBER(v_ice(:,:))
179	ENDIF
180	CALL lim_trp( kt ) ! -- Ice transport (Advection/diffusion)
181	IF( nn_limdyn == 2 .AND. nn_monocat /= 2 ) & ! -- Mechanical redistribution (ridging/rafting)
182	& CALL lim_itd_me
183	IF( nn_limdyn == 2 ) CALL lim_update1( kt ) ! -- Corrections
184	!
185	ENDIF
186
187	! ---
188	#if defined key_agrif
189	IF( .NOT. Agrif_Root() ) CALL agrif_interp_lim3('T')
190	#endif
191	IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice_lim( kt ) ! -- bdy ice thermo
192	! previous lead fraction and ice volume for flux calculations
193	CALL lim_var_glo2eqv ! ht_i and ht_s for ice albedo calculation
194	CALL lim_var_agg(1) ! at_i for coupling
195	CALL ice_bef
196
197	!------------------------------------------------------!
198	! --- Thermodynamical coupling with the atmosphere --- !
199	!------------------------------------------------------!
200	! It provides the following fields:
201	! qsr_ice , qns_ice : solar & non solar heat flux over ice (T-point) [W/m2]
202	! qla_ice : latent heat flux over ice (T-point) [W/m2]
203	! dqns_ice, dqla_ice : non solar & latent heat sensistivity (T-point) [W/m2]
204	! tprecip , sprecip : total & solid precipitation (T-point) [Kg/m2/s]
205	! fr1_i0 , fr2_i0 : 1sr & 2nd fraction of qsr penetration in ice [%]
206	!----------------------------------------------------------------------------------------
207
208	CALL albedo_ice( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos MV MP 2016
209
210	SELECT CASE( ksbc )
211	CASE( jp_usr ) ; CALL usrdef_sbc_ice_flx( kt ) ! user defined formulation
212	CASE( jp_blk ) ! bulk formulation
213	! albedo depends on cloud fraction because of non-linear spectral effects
214	alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
215	CALL blk_ice_flx( t_su, alb_ice )
216	IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
217	IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
218	CASE ( jp_purecpl )
219	! albedo depends on cloud fraction because of non-linear spectral effects
220	alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
221	CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
222	IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
223	END SELECT
224
225	!----------------------------!
226	! --- ice thermodynamics --- !
227	!----------------------------!
228	! --- zap this if no ice thermo --- !
229	IF( ln_limthd ) CALL lim_thd( kt ) ! -- Ice thermodynamics
230
231	! MV MP 2016
232	IF ( ln_pnd ) CALL lim_mp( kt ) ! -- Melt ponds
233	! END MV MP 2016
234
235	IF( ln_limthd ) CALL lim_update2( kt ) ! -- Corrections
236	! ---
237	# if defined key_agrif
238	IF( .NOT. Agrif_Root() ) CALL agrif_update_lim3( kt )
239	# endif
240	CALL lim_var_glo2eqv ! necessary calls (at least for coupling)
241	CALL lim_var_agg( 2 ) ! necessary calls (at least for coupling)
242	!
243	!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
244	!! moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
245	!!# if defined key_agrif
246	!! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid()
247	!!# endif
248	CALL lim_sbc_flx( kt ) ! -- Update surface ocean mass, heat and salt fluxes
249	!!# if defined key_agrif
250	!! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid()
251	!!# endif
252	IF( ln_limdiahsb ) CALL lim_diahsb( kt ) ! -- Diagnostics and outputs
253	!
254	CALL lim_wri( 1 ) ! -- Ice outputs
255	!
256	IF( kt == nit000 .AND. ln_rstart ) &
257	& CALL iom_close( numrir ) ! close input ice restart file
258	!
259	IF( lrst_ice ) CALL lim_rst_write( kt ) ! -- Ice restart file
260	!
261	IF( ln_limctl ) CALL lim_ctl( kt ) ! alerts in case of model crash
262	!
263	ENDIF ! End sea-ice time step only
264
265	!-------------------------!
266	! --- Ocean time step --- !
267	!-------------------------!
268	! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
269	! using before instantaneous surf. currents
270	IF( ln_limdyn ) CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )
271	!!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!!
272	!
273	IF( nn_timing == 1 ) CALL timing_stop('ice_stp')
274	!
275	END SUBROUTINE ice_stp
276
277
278	SUBROUTINE ice_init
279	!!----------------------------------------------------------------------
280	!! * ROUTINE ice_init *
281	!!
282	!! ** purpose : Allocate all the dynamic arrays of the LIM-3 modules
283	!!----------------------------------------------------------------------
284	INTEGER :: ji, jj, ierr
285	!!----------------------------------------------------------------------
286	IF(lwp) WRITE(numout,*)
287	IF(lwp) WRITE(numout,*) 'ice_init : update ocean surface boudary condition'
288	IF(lwp) WRITE(numout,*) '~~~~~~~~ via Louvain la Neuve Ice Model (LIM-3) time stepping'
289	!
290	! ! Open the reference and configuration namelist files and namelist output file
291	CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
292	CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
293	IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
294	!
295	CALL ice_run_init ! set some ice run parameters
296	!
297	! ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
298	ierr = ice_alloc () ! ice variables
299	ierr = ierr + sbc_ice_alloc () ! surface forcing
300	ierr = ierr + thd_ice_alloc () ! thermodynamics
301	IF( ln_limdyn ) ierr = ierr + lim_itd_me_alloc () ! ice thickness distribution - mechanics
302	!
303	IF( lk_mpp ) CALL mpp_sum( ierr )
304	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
305	!
306	CALL lim_dyn_init ! set ice dynamics parameters
307	!
308	CALL ice_itd_init ! ice thickness distribution initialization
309	!
310	CALL lim_thd_init ! set ice thermodynics parameters
311	!
312	CALL lim_thd_sal_init ! set ice salinity parameters
313
314	! MV MP 2016
315	CALL lim_mp_init ! set melt ponds parameters
316	! END MV MP 2016
317
318	IF( ln_limdyn ) CALL lim_itd_me_init ! ice thickness distribution initialization for mecanical deformation
319	! ! Initial sea-ice state
320	IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst
321	CALL lim_istate
322	ELSE ! start from a restart file
323	CALL lim_rst_read
324	ENDIF
325	CALL lim_var_agg(2)
326	CALL lim_var_glo2eqv
327	!
328	CALL lim_sbc_init ! ice surface boundary condition
329	!
330	IF( ln_limdiahsb) CALL lim_diahsb_init ! initialization for diags
331	!
332	fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction
333	tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu
334	!
335	DO jj = 1, jpj
336	DO ji = 1, jpi
337	IF( gphit(ji,jj) > 0._wp ) THEN ; rn_amax_2d(ji,jj) = rn_amax_n ! NH
338	ELSE ; rn_amax_2d(ji,jj) = rn_amax_s ! SH
339	ENDIF
340	END DO
341	END DO
342	!
343	END SUBROUTINE ice_init
344
345
346	SUBROUTINE ice_run_init
347	!!-------------------------------------------------------------------
348	!! * ROUTINE ice_run_init *
349	!!
350	!! ** Purpose : Definition some run parameter for ice model
351	!!
352	!! ** Method : Read the namicerun namelist and check the parameter
353	!! values called at the first timestep (nit000)
354	!!
355	!! ** input : Namelist namicerun
356	!!-------------------------------------------------------------------
357	INTEGER :: ios ! Local integer output status for namelist read
358	NAMELIST/namicerun/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir, &
359	& cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice
360	NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt
361	!!-------------------------------------------------------------------
362	!
363	REWIND( numnam_ice_ref ) ! Namelist namicerun in reference namelist : Parameters for ice
364	READ ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
365	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
366
367	REWIND( numnam_ice_cfg ) ! Namelist namicerun in configuration namelist : Parameters for ice
368	READ ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
369	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
370	IF(lwm) WRITE ( numoni, namicerun )
371	!
372	REWIND( numnam_ice_ref ) ! Namelist namicediag in reference namelist : Parameters for ice
373	READ ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
374	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
375
376	REWIND( numnam_ice_cfg ) ! Namelist namicediag in configuration namelist : Parameters for ice
377	READ ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
378	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
379	IF(lwm) WRITE ( numoni, namicediag )
380	!
381	IF(lwp) THEN ! control print
382	WRITE(numout,*)
383	WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
384	WRITE(numout,*) ' ~~~~~~'
385	WRITE(numout,*) ' number of ice categories jpl = ', jpl
386	WRITE(numout,*) ' number of ice layers nlay_i = ', nlay_i
387	WRITE(numout,*) ' number of snow layers nlay_s = ', nlay_s
388	WRITE(numout,*) ' virtual ITD mono-category param (1-4) or not (0) nn_monocat = ', nn_monocat
389	WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n
390	WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s
391	WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_limthd = ', ln_limthd
392	WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_limdyn = ', ln_limdyn
393	WRITE(numout,*) ' (ln_limdyn=T) Ice dynamics switch nn_limdyn = ', nn_limdyn
394	WRITE(numout,*) ' 2: total'
395	WRITE(numout,*) ' 1: advection only (no diffusion, no ridging/rafting)'
396	WRITE(numout,*) ' 0: advection only (as 1 + prescribed velocity, bypass rheology)'
397	WRITE(numout,*) ' (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0) = ', rn_uice
398	WRITE(numout,*) ' (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0) = ', rn_vice
399	WRITE(numout,*)
400	WRITE(numout,*) '...and ice diagnostics'
401	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
402	WRITE(numout,*) ' Diagnose online heat/mass/salt budget ln_limdiachk = ', ln_limdiachk
403	WRITE(numout,*) ' Output heat/mass/salt budget ln_limdiahsb = ', ln_limdiahsb
404	WRITE(numout,*) ' control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl
405	WRITE(numout,*) ' i-index for control prints (ln_limctl=true) = ', iiceprt
406	WRITE(numout,*) ' j-index for control prints (ln_limctl=true) = ', jiceprt
407	ENDIF
408	!
409	IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
410	nn_monocat = 0
411	IF(lwp) WRITE(numout,*)
412	IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
413	ENDIF
414	IF ( ( jpl == 1 ) .AND. ( nn_monocat == 0 ) ) THEN
415	CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
416	ENDIF
417	!
418	! sea-ice timestep and inverse
419	rdt_ice = REAL(nn_fsbc) * rdt
420	r1_rdtice = 1._wp / rdt_ice
421
422	! inverse of nlay_i and nlay_s
423	r1_nlay_i = 1._wp / REAL( nlay_i, wp )
424	r1_nlay_s = 1._wp / REAL( nlay_s, wp )
425	!
426	IF( lwp .AND. ln_bdy .AND. ln_limdiachk ) &
427	& CALL ctl_warn('online conservation check activated but it does not work with BDY')
428	!
429	IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice
430	!
431	END SUBROUTINE ice_run_init
432
433
434	SUBROUTINE ice_itd_init
435	!!------------------------------------------------------------------
436	!! * ROUTINE ice_itd_init *
437	!!
438	!! ** Purpose : Initializes the ice thickness distribution
439	!! ** Method : ...
440	!! ** input : Namelist namiceitd
441	!!-------------------------------------------------------------------
442	INTEGER :: ios ! Local integer output status for namelist read
443	NAMELIST/namiceitd/ rn_himean
444	!
445	INTEGER :: jl ! dummy loop index
446	REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars
447	REAL(wp) :: zhmax, znum, zden, zalpha !
448	!!------------------------------------------------------------------
449	!
450	REWIND( numnam_ice_ref ) ! Namelist namiceitd in reference namelist : Parameters for ice
451	READ ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
452	905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
453
454	REWIND( numnam_ice_cfg ) ! Namelist namiceitd in configuration namelist : Parameters for ice
455	READ ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
456	906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
457	IF(lwm) WRITE ( numoni, namiceitd )
458	!
459	IF(lwp) THEN ! control print
460	WRITE(numout,*)
461	WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution '
462	WRITE(numout,*) '~~~~~~~~~~~~'
463	WRITE(numout,*) ' mean ice thickness in the domain rn_himean = ', rn_himean
464	ENDIF
465	!
466	!----------------------------------
467	!- Thickness categories boundaries
468	!----------------------------------
469	!
470	!== h^(-alpha) function ==!
471	zalpha = 0.05_wp
472	zhmax = 3._wp * rn_himean
473	DO jl = 1, jpl
474	znum = jpl * ( zhmax+1 )**zalpha
475	zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
476	hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
477	END DO
478	!
479	DO jl = 1, jpl ! mean thickness by category
480	hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
481	END DO
482	!
483	hi_max(jpl) = 99._wp ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
484	!
485	IF(lwp) WRITE(numout,*)
486	IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
487	IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
488	!
489	END SUBROUTINE ice_itd_init
490
491
492	SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
493	!!---------------------------------------------------------------------
494	!! * ROUTINE ice_lim_flx *
495	!!
496	!! ** Purpose : update the ice surface boundary condition by averaging and / or
497	!! redistributing fluxes on ice categories
498	!!
499	!! ** Method : average then redistribute
500	!!
501	!! ** Action :
502	!!---------------------------------------------------------------------
503	INTEGER , INTENT(in ) :: k_limflx ! =-1 do nothing; =0 average ;
504	! ! = 1 average and redistribute ; =2 redistribute
505	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature
506	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo
507	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqns_ice ! non solar flux
508	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux
509	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity
510	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation
511	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity
512	!
513	INTEGER :: jl ! dummy loop index
514	!
515	REAL(wp), DIMENSION(jpi,jpj) :: zalb_m ! Mean albedo over all categories
516	REAL(wp), DIMENSION(jpi,jpj) :: ztem_m ! Mean temperature over all categories
517	!
518	REAL(wp), DIMENSION(jpi,jpj) :: z_qsr_m ! Mean solar heat flux over all categories
519	REAL(wp), DIMENSION(jpi,jpj) :: z_qns_m ! Mean non solar heat flux over all categories
520	REAL(wp), DIMENSION(jpi,jpj) :: z_evap_m ! Mean sublimation over all categories
521	REAL(wp), DIMENSION(jpi,jpj) :: z_dqn_m ! Mean d(qns)/dT over all categories
522	REAL(wp), DIMENSION(jpi,jpj) :: z_devap_m ! Mean d(evap)/dT over all categories
523	!!----------------------------------------------------------------------
524	!
525	IF( nn_timing == 1 ) CALL timing_start('ice_lim_flx')
526	!
527	SELECT CASE( k_limflx ) !== averaged on all ice categories ==!
528	CASE( 0 , 1 )
529	!
530	z_qns_m (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
531	z_qsr_m (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
532	z_dqn_m (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
533	z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
534	z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
535	DO jl = 1, jpl
536	pdqn_ice (:,:,jl) = z_dqn_m(:,:)
537	pdevap_ice(:,:,jl) = z_devap_m(:,:)
538	END DO
539	!
540	DO jl = 1, jpl
541	pqns_ice (:,:,jl) = z_qns_m(:,:)
542	pqsr_ice (:,:,jl) = z_qsr_m(:,:)
543	pevap_ice(:,:,jl) = z_evap_m(:,:)
544	END DO
545	!
546	END SELECT
547	!
548	SELECT CASE( k_limflx ) !== redistribution on all ice categories ==!
549	CASE( 1 , 2 )
550	!
551	zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )
552	ztem_m(:,:) = fice_ice_ave ( ptn_ice (:,:,:) )
553	DO jl = 1, jpl
554	pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
555	pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
556	pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
557	END DO
558	!
559	END SELECT
560	!
561	IF( nn_timing == 1 ) CALL timing_stop('ice_lim_flx')
562	!
563	END SUBROUTINE ice_lim_flx
564
565
566	SUBROUTINE ice_bef
567	!!----------------------------------------------------------------------
568	!! * ROUTINE ice_bef *
569	!!
570	!! ** purpose : store ice variables at "before" time step
571	!!----------------------------------------------------------------------
572	INTEGER :: ji, jj, jl ! dummy loop index
573	!!----------------------------------------------------------------------
574	!
575	DO jl = 1, jpl
576
577	DO jj = 1, jpj
578	DO ji = 1, jpi
579	a_i_b (ji,jj,jl) = a_i (ji,jj,jl) ! ice area
580	v_i_b (ji,jj,jl) = v_i (ji,jj,jl) ! ice volume
581	v_s_b (ji,jj,jl) = v_s (ji,jj,jl) ! snow volume
582	smv_i_b(ji,jj,jl) = smv_i(ji,jj,jl) ! salt content
583	oa_i_b (ji,jj,jl) = oa_i (ji,jj,jl) ! areal age content
584	e_s_b (ji,jj,:,jl) = e_s (ji,jj,:,jl) ! snow thermal energy
585	e_i_b (ji,jj,:,jl) = e_i (ji,jj,:,jl) ! ice thermal energy
586	! ! ice thickness
587	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes
588	ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
589	ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
590	END DO
591	END DO
592
593	END DO
594
595	! ice velocities & total concentration
596	DO jj = 1, jpj
597	DO ji = 1, jpi
598	at_i_b(ji,jj) = SUM( a_i_b(ji,jj,:) )
599	u_ice_b(ji,jj) = u_ice(ji,jj)
600	v_ice_b(ji,jj) = v_ice(ji,jj)
601	END DO
602	END DO
603
604	END SUBROUTINE ice_bef
605
606
607	SUBROUTINE ice_diag0
608	!!----------------------------------------------------------------------
609	!! * ROUTINE ice_diag0 *
610	!!
611	!! ** purpose : set ice-ocean and ice-atm. fluxes to zeros at the beggining
612	!! of the time step
613	!!----------------------------------------------------------------------
614	INTEGER :: ji, jj ! dummy loop index
615	!!----------------------------------------------------------------------
616	DO jj = 1, jpj
617	DO ji = 1, jpi
618	sfx (ji,jj) = 0._wp ;
619	sfx_bri(ji,jj) = 0._wp ; sfx_lam(ji,jj) = 0._wp
620	sfx_sni(ji,jj) = 0._wp ; sfx_opw(ji,jj) = 0._wp
621	sfx_bog(ji,jj) = 0._wp ; sfx_dyn(ji,jj) = 0._wp
622	sfx_bom(ji,jj) = 0._wp ; sfx_sum(ji,jj) = 0._wp
623	sfx_res(ji,jj) = 0._wp ; sfx_sub(ji,jj) = 0._wp
624	!
625	wfx_snw(ji,jj) = 0._wp ; wfx_ice(ji,jj) = 0._wp
626	wfx_sni(ji,jj) = 0._wp ; wfx_opw(ji,jj) = 0._wp
627	wfx_bog(ji,jj) = 0._wp ; wfx_dyn(ji,jj) = 0._wp
628	wfx_bom(ji,jj) = 0._wp ; wfx_sum(ji,jj) = 0._wp
629	wfx_res(ji,jj) = 0._wp ; wfx_sub(ji,jj) = 0._wp
630	wfx_spr(ji,jj) = 0._wp ; wfx_lam(ji,jj) = 0._wp
631	wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
632	wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
633	wfx_snw_sni(ji,jj) = 0._wp
634	! MV MP 2016
635	wfx_pnd(ji,jj) = 0._wp
636	! END MV MP 2016
637
638	hfx_thd(ji,jj) = 0._wp ;
639	hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp
640	hfx_bog(ji,jj) = 0._wp ; hfx_dyn(ji,jj) = 0._wp
641	hfx_bom(ji,jj) = 0._wp ; hfx_sum(ji,jj) = 0._wp
642	hfx_res(ji,jj) = 0._wp ; hfx_sub(ji,jj) = 0._wp
643	hfx_spr(ji,jj) = 0._wp ; hfx_dif(ji,jj) = 0._wp
644	hfx_err(ji,jj) = 0._wp ; hfx_err_rem(ji,jj) = 0._wp
645	hfx_err_dif(ji,jj) = 0._wp
646	wfx_err_sub(ji,jj) = 0._wp
647	!
648	afx_tot(ji,jj) = 0._wp ;
649	afx_dyn(ji,jj) = 0._wp ; afx_thd(ji,jj) = 0._wp
650	!
651	diag_heat(ji,jj) = 0._wp ; diag_smvi(ji,jj) = 0._wp
652	diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp
653
654	! SIMIP diagnostics
655	diag_dms_dyn(ji,jj) = 0._wp ; diag_dmi_dyn(ji,jj) = 0._wp
656	diag_fc_bo(ji,jj) = 0._wp ; diag_fc_su(ji,jj) = 0._wp
657
658	tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
659	END DO
660	END DO
661
662	END SUBROUTINE ice_diag0
663
664
665	FUNCTION fice_cell_ave ( ptab )
666	!!--------------------------------------------------------------------------
667	!! * Compute average over categories, for grid cell (ice covered and free ocean)
668	!!--------------------------------------------------------------------------
669	REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
670	REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
671	INTEGER :: jl ! Dummy loop index
672
673	fice_cell_ave (:,:) = 0._wp
674	DO jl = 1, jpl
675	fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
676	END DO
677
678	END FUNCTION fice_cell_ave
679
680
681	FUNCTION fice_ice_ave ( ptab )
682	!!--------------------------------------------------------------------------
683	!! * Compute average over categories, for ice covered part of grid cell
684	!!--------------------------------------------------------------------------
685	REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
686	REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
687
688	fice_ice_ave (:,:) = 0.0_wp
689	WHERE ( at_i (:,:) > 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
690
691	END FUNCTION fice_ice_ave
692
693	#else
694	!!----------------------------------------------------------------------
695	!! Default option Dummy module NO LIM 3.0 sea-ice model
696	!!----------------------------------------------------------------------
697	CONTAINS
698	!
699	SUBROUTINE ice_stp ( kt, ksbc ) ! Dummy routine
700	INTEGER, INTENT(in) :: kt, ksbc
701	WRITE(,) 'ice_stp: You should not have seen this print! error?', kt
702	END SUBROUTINE ice_stp
703	!
704	SUBROUTINE ice_init ! Dummy routine
705	END SUBROUTINE ice_init
706	#endif
707
708	!!======================================================================
709	END MODULE icestp

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

Download in other formats: