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

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

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

Download in other formats: