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

Last change on this file since 8341 was 8341, checked in by clem, 7 years ago
correct the heat conservation (all fine except limthd_da for a reason I do not understand)
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 + sbc_ice_alloc () ! surface forcing
307	ierr = ierr + thd_ice_alloc () ! thermodynamics
308	IF( ln_limdyn ) ierr = ierr + lim_itd_me_alloc () ! ice thickness distribution - mechanics
309	!
310	IF( lk_mpp ) CALL mpp_sum( ierr )
311	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
312	!
313	CALL lim_dyn_init ! set ice dynamics parameters
314	!
315	CALL ice_itd_init ! ice thickness distribution initialization
316	!
317	CALL lim_thd_init ! set ice thermodynics parameters
318	!
319	CALL lim_thd_sal_init ! set ice salinity parameters
320
321	! MV MP 2016
322	CALL lim_mp_init ! set melt ponds parameters
323	! END MV MP 2016
324
325	IF( ln_limdyn ) CALL lim_itd_me_init ! ice thickness distribution initialization for mecanical deformation
326	! ! Initial sea-ice state
327	IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst
328	CALL lim_istate
329	ELSE ! start from a restart file
330	CALL lim_rst_read
331	ENDIF
332	CALL lim_var_agg(2)
333	CALL lim_var_glo2eqv
334	!
335	CALL lim_sbc_init ! ice surface boundary condition
336	!
337	IF( ln_limdiahsb) CALL lim_diahsb_init ! initialization for diags
338	!
339	fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction
340	tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu
341	!
342	DO jj = 1, jpj
343	DO ji = 1, jpi
344	IF( gphit(ji,jj) > 0._wp ) THEN ; rn_amax_2d(ji,jj) = rn_amax_n ! NH
345	ELSE ; rn_amax_2d(ji,jj) = rn_amax_s ! SH
346	ENDIF
347	END DO
348	END DO
349	!
350	END SUBROUTINE ice_init
351
352
353	SUBROUTINE ice_run_init
354	!!-------------------------------------------------------------------
355	!! * ROUTINE ice_run_init *
356	!!
357	!! ** Purpose : Definition some run parameter for ice model
358	!!
359	!! ** Method : Read the namicerun namelist and check the parameter
360	!! values called at the first timestep (nit000)
361	!!
362	!! ** input : Namelist namicerun
363	!!-------------------------------------------------------------------
364	INTEGER :: ios ! Local integer output status for namelist read
365	NAMELIST/namicerun/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir, &
366	& cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice
367	NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt
368	!!-------------------------------------------------------------------
369	!
370	REWIND( numnam_ice_ref ) ! Namelist namicerun in reference namelist : Parameters for ice
371	READ ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
372	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
373
374	REWIND( numnam_ice_cfg ) ! Namelist namicerun in configuration namelist : Parameters for ice
375	READ ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
376	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
377	IF(lwm) WRITE ( numoni, namicerun )
378	!
379	REWIND( numnam_ice_ref ) ! Namelist namicediag in reference namelist : Parameters for ice
380	READ ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
381	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
382
383	REWIND( numnam_ice_cfg ) ! Namelist namicediag in configuration namelist : Parameters for ice
384	READ ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
385	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
386	IF(lwm) WRITE ( numoni, namicediag )
387	!
388	IF(lwp) THEN ! control print
389	WRITE(numout,*)
390	WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
391	WRITE(numout,*) ' ~~~~~~'
392	WRITE(numout,*) ' number of ice categories jpl = ', jpl
393	WRITE(numout,*) ' number of ice layers nlay_i = ', nlay_i
394	WRITE(numout,*) ' number of snow layers nlay_s = ', nlay_s
395	WRITE(numout,*) ' virtual ITD mono-category param (1-4) or not (0) nn_monocat = ', nn_monocat
396	WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n
397	WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s
398	WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_limthd = ', ln_limthd
399	WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_limdyn = ', ln_limdyn
400	WRITE(numout,*) ' (ln_limdyn=T) Ice dynamics switch nn_limdyn = ', nn_limdyn
401	WRITE(numout,*) ' 2: total'
402	WRITE(numout,*) ' 1: advection only (no diffusion, no ridging/rafting)'
403	WRITE(numout,*) ' 0: advection only (as 1 + prescribed velocity, bypass rheology)'
404	WRITE(numout,*) ' (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0) = ', rn_uice
405	WRITE(numout,*) ' (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0) = ', rn_vice
406	WRITE(numout,*)
407	WRITE(numout,*) '...and ice diagnostics'
408	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
409	WRITE(numout,*) ' Diagnose online heat/mass/salt budget ln_limdiachk = ', ln_limdiachk
410	WRITE(numout,*) ' Output heat/mass/salt budget ln_limdiahsb = ', ln_limdiahsb
411	WRITE(numout,*) ' control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl
412	WRITE(numout,*) ' i-index for control prints (ln_limctl=true) = ', iiceprt
413	WRITE(numout,*) ' j-index for control prints (ln_limctl=true) = ', jiceprt
414	ENDIF
415	!
416	IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
417	nn_monocat = 0
418	IF(lwp) WRITE(numout,*)
419	IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
420	ENDIF
421	IF ( ( jpl == 1 ) .AND. ( nn_monocat == 0 ) ) THEN
422	CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
423	ENDIF
424	!
425	! sea-ice timestep and inverse
426	rdt_ice = REAL(nn_fsbc) * rdt
427	r1_rdtice = 1._wp / rdt_ice
428
429	! inverse of nlay_i and nlay_s
430	r1_nlay_i = 1._wp / REAL( nlay_i, wp )
431	r1_nlay_s = 1._wp / REAL( nlay_s, wp )
432	!
433	IF( lwp .AND. ln_bdy .AND. ln_limdiachk ) &
434	& CALL ctl_warn('online conservation check activated but it does not work with BDY')
435	!
436	IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice
437	!
438	END SUBROUTINE ice_run_init
439
440
441	SUBROUTINE ice_itd_init
442	!!------------------------------------------------------------------
443	!! * ROUTINE ice_itd_init *
444	!!
445	!! ** Purpose : Initializes the ice thickness distribution
446	!! ** Method : ...
447	!! ** input : Namelist namiceitd
448	!!-------------------------------------------------------------------
449	INTEGER :: ios ! Local integer output status for namelist read
450	NAMELIST/namiceitd/ rn_himean
451	!
452	INTEGER :: jl ! dummy loop index
453	REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars
454	REAL(wp) :: zhmax, znum, zden, zalpha !
455	!!------------------------------------------------------------------
456	!
457	REWIND( numnam_ice_ref ) ! Namelist namiceitd in reference namelist : Parameters for ice
458	READ ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
459	905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
460
461	REWIND( numnam_ice_cfg ) ! Namelist namiceitd in configuration namelist : Parameters for ice
462	READ ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
463	906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
464	IF(lwm) WRITE ( numoni, namiceitd )
465	!
466	IF(lwp) THEN ! control print
467	WRITE(numout,*)
468	WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution '
469	WRITE(numout,*) '~~~~~~~~~~~~'
470	WRITE(numout,*) ' mean ice thickness in the domain rn_himean = ', rn_himean
471	ENDIF
472	!
473	!----------------------------------
474	!- Thickness categories boundaries
475	!----------------------------------
476	!
477	hi_max(:) = 0._wp
478	!
479	!== h^(-alpha) function ==!
480	zalpha = 0.05_wp
481	zhmax = 3._wp * rn_himean
482	DO jl = 1, jpl
483	znum = jpl * ( zhmax+1 )**zalpha
484	zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
485	hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
486	END DO
487	!
488	!
489	DO jl = 1, jpl ! mean thickness by category
490	hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
491	END DO
492	!
493	hi_max(jpl) = 99._wp ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
494	!
495	IF(lwp) WRITE(numout,*)
496	IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
497	IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
498	!
499	END SUBROUTINE ice_itd_init
500
501
502	SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
503	!!---------------------------------------------------------------------
504	!! * ROUTINE ice_lim_flx *
505	!!
506	!! ** Purpose : update the ice surface boundary condition by averaging and / or
507	!! redistributing fluxes on ice categories
508	!!
509	!! ** Method : average then redistribute
510	!!
511	!! ** Action :
512	!!---------------------------------------------------------------------
513	INTEGER , INTENT(in ) :: k_limflx ! =-1 do nothing; =0 average ;
514	! ! = 1 average and redistribute ; =2 redistribute
515	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature
516	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo
517	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqns_ice ! non solar flux
518	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux
519	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity
520	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation
521	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity
522	!
523	INTEGER :: jl ! dummy loop index
524	!
525	REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m ! Mean albedo over all categories
526	REAL(wp), POINTER, DIMENSION(:,:) :: ztem_m ! Mean temperature over all categories
527	!
528	REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories
529	REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m ! Mean non solar heat flux over all categories
530	REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m ! Mean sublimation over all categories
531	REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m ! Mean d(qns)/dT over all categories
532	REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories
533	!!----------------------------------------------------------------------
534	!
535	IF( nn_timing == 1 ) CALL timing_start('ice_lim_flx')
536	!
537	SELECT CASE( k_limflx ) !== averaged on all ice categories ==!
538	CASE( 0 , 1 )
539	CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
540	!
541	z_qns_m (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
542	z_qsr_m (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
543	z_dqn_m (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
544	z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
545	z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
546	DO jl = 1, jpl
547	pdqn_ice (:,:,jl) = z_dqn_m(:,:)
548	pdevap_ice(:,:,jl) = z_devap_m(:,:)
549	END DO
550	!
551	DO jl = 1, jpl
552	pqns_ice (:,:,jl) = z_qns_m(:,:)
553	pqsr_ice (:,:,jl) = z_qsr_m(:,:)
554	pevap_ice(:,:,jl) = z_evap_m(:,:)
555	END DO
556	!
557	CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
558	END SELECT
559	!
560	SELECT CASE( k_limflx ) !== redistribution on all ice categories ==!
561	CASE( 1 , 2 )
562	CALL wrk_alloc( jpi,jpj, zalb_m, ztem_m )
563	!
564	zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )
565	ztem_m(:,:) = fice_ice_ave ( ptn_ice (:,:,:) )
566	DO jl = 1, jpl
567	pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
568	pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
569	pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
570	END DO
571	!
572	CALL wrk_dealloc( jpi,jpj, zalb_m, ztem_m )
573	END SELECT
574	!
575	IF( nn_timing == 1 ) CALL timing_stop('ice_lim_flx')
576	!
577	END SUBROUTINE ice_lim_flx
578
579
580	SUBROUTINE ice_bef
581	!!----------------------------------------------------------------------
582	!! * ROUTINE ice_bef *
583	!!
584	!! ** purpose : store ice variables at "before" time step
585	!!----------------------------------------------------------------------
586	a_i_b (:,:,:) = a_i (:,:,:) ! ice area
587	e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy
588	v_i_b (:,:,:) = v_i (:,:,:) ! ice volume
589	v_s_b (:,:,:) = v_s (:,:,:) ! snow volume
590	e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy
591	smv_i_b(:,:,:) = smv_i(:,:,:) ! salt content
592	oa_i_b (:,:,:) = oa_i (:,:,:) ! areal age content
593	u_ice_b(:,:) = u_ice(:,:)
594	v_ice_b(:,:) = v_ice(:,:)
595	!
596	at_i_b (:,:) = SUM( a_i_b(:,:,:), dim=3 )
597
598	END SUBROUTINE ice_bef
599
600
601	SUBROUTINE ice_diag0
602	!!----------------------------------------------------------------------
603	!! * ROUTINE ice_diag0 *
604	!!
605	!! ** purpose : set ice-ocean and ice-atm. fluxes to zeros at the beggining
606	!! of the time step
607	!!----------------------------------------------------------------------
608	sfx (:,:) = 0._wp ;
609	sfx_bri(:,:) = 0._wp ; sfx_lam(:,:) = 0._wp
610	sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp
611	sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp
612	sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp
613	sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp
614	!
615	wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp
616	wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp
617	wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp
618	wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp
619	wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp
620	wfx_spr(:,:) = 0._wp ; wfx_lam(:,:) = 0._wp
621	wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
622	wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
623	wfx_snw_sni(:,:) = 0._wp
624	! MV MP 2016
625	wfx_pnd(:,:) = 0._wp
626	! END MV MP 2016
627
628	hfx_thd(:,:) = 0._wp ;
629	hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp
630	hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp
631	hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp
632	hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp
633	hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp
634	hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp
635	hfx_err_dif(:,:) = 0._wp
636	wfx_err_sub(:,:) = 0._wp
637	!
638	afx_tot(:,:) = 0._wp ;
639	afx_dyn(:,:) = 0._wp ; afx_thd(:,:) = 0._wp
640	!
641	diag_heat(:,:) = 0._wp ; diag_smvi(:,:) = 0._wp
642	diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp
643
644	! SIMIP diagnostics
645	diag_dms_dyn(:,:) = 0._wp ; diag_dmi_dyn(:,:) = 0._wp
646	diag_fc_bo(:,:) = 0._wp ; diag_fc_su(:,:) = 0._wp
647
648	tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
649
650	END SUBROUTINE ice_diag0
651
652
653	FUNCTION fice_cell_ave ( ptab )
654	!!--------------------------------------------------------------------------
655	!! * Compute average over categories, for grid cell (ice covered and free ocean)
656	!!--------------------------------------------------------------------------
657	REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
658	REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
659	INTEGER :: jl ! Dummy loop index
660
661	fice_cell_ave (:,:) = 0._wp
662	DO jl = 1, jpl
663	fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
664	END DO
665
666	END FUNCTION fice_cell_ave
667
668
669	FUNCTION fice_ice_ave ( ptab )
670	!!--------------------------------------------------------------------------
671	!! * Compute average over categories, for ice covered part of grid cell
672	!!--------------------------------------------------------------------------
673	REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
674	REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
675
676	fice_ice_ave (:,:) = 0.0_wp
677	WHERE ( at_i (:,:) > 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
678
679	END FUNCTION fice_ice_ave
680
681	#else
682	!!----------------------------------------------------------------------
683	!! Default option Dummy module NO LIM 3.0 sea-ice model
684	!!----------------------------------------------------------------------
685	CONTAINS
686	!
687	SUBROUTINE ice_stp ( kt, ksbc ) ! Dummy routine
688	INTEGER, INTENT(in) :: kt, ksbc
689	WRITE(,) 'ice_stp: You should not have seen this print! error?', kt
690	END SUBROUTINE ice_stp
691	!
692	SUBROUTINE ice_init ! Dummy routine
693	END SUBROUTINE ice_init
694	#endif
695
696	!!======================================================================
697	END MODULE icestp

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

Download in other formats: