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.

sbcice_lim.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 8316

Last change on this file since 8316 was 8316, checked in by clem, 7 years ago
STEP2 (3): remove obsolete features
Property svn:keywords set to `Id`
File size: 34.5 KB

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

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

Download in other formats: