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

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

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

Download in other formats: