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

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

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

Download in other formats: