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

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

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

Download in other formats: