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/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 5051

Last change on this file since 5051 was 5051, checked in by clem, 9 years ago
LIM3 initialization is now called at the same time as other sbc fields
Property svn:keywords set to `Id`
File size: 32.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	!!----------------------------------------------------------------------
16	#if defined key_lim3
17	!!----------------------------------------------------------------------
18	!! 'key_lim3' : LIM 3.0 sea-ice model
19	!!----------------------------------------------------------------------
20	!! sbc_ice_lim : sea-ice model time-stepping and update ocean sbc over ice-covered area
21	!!----------------------------------------------------------------------
22	USE oce ! ocean dynamics and tracers
23	USE dom_oce ! ocean space and time domain
24	USE par_ice ! sea-ice parameters
25	USE ice ! LIM-3: ice variables
26	USE thd_ice ! LIM-3: thermodynamical variables
27	USE dom_ice ! LIM-3: ice domain
28
29	USE sbc_oce ! Surface boundary condition: ocean fields
30	USE sbc_ice ! Surface boundary condition: ice fields
31	USE sbcblk_core ! Surface boundary condition: CORE bulk
32	USE sbcblk_clio ! Surface boundary condition: CLIO bulk
33	USE sbccpl ! Surface boundary condition: coupled interface
34	USE albedo ! ocean & ice albedo
35
36	USE phycst ! Define parameters for the routines
37	USE eosbn2 ! equation of state
38	USE limdyn ! Ice dynamics
39	USE limtrp ! Ice transport
40	USE limthd ! Ice thermodynamics
41	USE limitd_th ! Thermodynamics on ice thickness distribution
42	USE limitd_me ! Mechanics on ice thickness distribution
43	USE limsbc ! sea surface boundary condition
44	USE limdiahsb ! Ice budget diagnostics
45	USE limwri ! Ice outputs
46	USE limrst ! Ice restarts
47	USE limupdate1 ! update of global variables
48	USE limupdate2 ! update of global variables
49	USE limvar ! Ice variables switch
50
51	USE limmsh ! LIM mesh
52	USE limistate ! LIM initial state
53	USE limthd_sal ! LIM ice thermodynamics: salinity
54
55	USE c1d ! 1D vertical configuration
56	USE lbclnk ! lateral boundary condition - MPP link
57	USE lib_mpp ! MPP library
58	USE wrk_nemo ! work arrays
59	USE timing ! Timing
60	USE iom ! I/O manager library
61	USE in_out_manager ! I/O manager
62	USE prtctl ! Print control
63	USE lib_fortran !
64	USE limctl
65
66	#if defined key_bdy
67	USE bdyice_lim ! unstructured open boundary data (bdy_ice_lim routine)
68	#endif
69
70	IMPLICIT NONE
71	PRIVATE
72
73	PUBLIC sbc_ice_lim ! routine called by sbcmod.F90
74	PUBLIC sbc_lim_init ! routine called by sbcmod.F90
75
76	!! * Substitutions
77	# include "domzgr_substitute.h90"
78	# include "vectopt_loop_substitute.h90"
79	!!----------------------------------------------------------------------
80	!! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
81	!! $Id$
82	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
83	!!----------------------------------------------------------------------
84	CONTAINS
85
86	!!======================================================================
87
88	SUBROUTINE sbc_ice_lim( kt, kblk )
89	!!---------------------------------------------------------------------
90	!! * ROUTINE sbc_ice_lim *
91	!!
92	!! ** Purpose : update the ocean surface boundary condition via the
93	!! Louvain la Neuve Sea Ice Model time stepping
94	!!
95	!! ** Method : ice model time stepping
96	!! - call the ice dynamics routine
97	!! - call the ice advection/diffusion routine
98	!! - call the ice thermodynamics routine
99	!! - call the routine that computes mass and
100	!! heat fluxes at the ice/ocean interface
101	!! - save the outputs
102	!! - save the outputs for restart when necessary
103	!!
104	!! ** Action : - time evolution of the LIM sea-ice model
105	!! - update all sbc variables below sea-ice:
106	!! utau, vtau, taum, wndm, qns , qsr, emp , sfx
107	!!---------------------------------------------------------------------
108	INTEGER, INTENT(in) :: kt ! ocean time step
109	INTEGER, INTENT(in) :: kblk ! type of bulk (=3 CLIO, =4 CORE, =5 COUPLED)
110	!!
111	INTEGER :: jl ! dummy loop index
112	REAL(wp) :: zcoef ! local scalar
113	REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_os, zalb_cs ! ice albedo under overcast/clear sky
114	REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_ice ! mean ice albedo (for coupled)
115	!!----------------------------------------------------------------------
116
117	IF( nn_timing == 1 ) CALL timing_start('sbc_ice_lim')
118
119	! !----------------------!
120	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Ice time-step only !
121	! !----------------------!
122	! ! Bulk Formulae !
123	! !----------------!
124	!
125	u_oce(:,:) = ssu_m(:,:) * umask(:,:,1) ! mean surface ocean current at ice velocity point
126	v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1) ! (C-grid dynamics : U- & V-points as the ocean)
127
128	! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
129	t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
130	!
131	! Ice albedo
132	CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
133	CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
134
135	SELECT CASE( kblk )
136	CASE( jp_core , jp_cpl ) ! CORE and COUPLED bulk formulations
137
138	! albedo depends on cloud fraction because of non-linear spectral effects
139	zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
140	! In CLIO the cloud fraction is read in the climatology and the all-sky albedo
141	! (zalb_ice) is computed within the bulk routine
142
143	END SELECT
144
145	! Mask sea ice surface temperature (set to rt0 over land)
146	DO jl = 1, jpl
147	t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
148	END DO
149
150	! Bulk formulae - provides the following fields:
151	! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
152	! qsr_ice , qns_ice : solar & non solar heat flux over ice (T-point) [W/m2]
153	! qla_ice : latent heat flux over ice (T-point) [W/m2]
154	! dqns_ice, dqla_ice : non solar & latent heat sensistivity (T-point) [W/m2]
155	! tprecip , sprecip : total & solid precipitation (T-point) [Kg/m2/s]
156	! fr1_i0 , fr2_i0 : 1sr & 2nd fraction of qsr penetration in ice [%]
157	!
158	SELECT CASE( kblk )
159	CASE( jp_clio ) ! CLIO bulk formulation
160	CALL blk_ice_clio( t_su , zalb_cs , zalb_os , zalb_ice , &
161	& utau_ice , vtau_ice , qns_ice , qsr_ice , &
162	& qla_ice , dqns_ice , dqla_ice , &
163	& tprecip , sprecip , &
164	& fr1_i0 , fr2_i0 , cp_ice_msh, jpl )
165	!
166	IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice , &
167	& dqns_ice, qla_ice, dqla_ice, nn_limflx )
168
169	CASE( jp_core ) ! CORE bulk formulation
170	CALL blk_ice_core( t_su , u_ice , v_ice , zalb_ice , &
171	& utau_ice , vtau_ice , qns_ice , qsr_ice , &
172	& qla_ice , dqns_ice , dqla_ice , &
173	& tprecip , sprecip , &
174	& fr1_i0 , fr2_i0 , cp_ice_msh, jpl )
175	!
176	IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice , &
177	& dqns_ice, qla_ice, dqla_ice, nn_limflx )
178	!
179	CASE ( jp_cpl )
180
181	CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
182
183	END SELECT
184
185	! !----------------------!
186	! ! LIM-3 time-stepping !
187	! !----------------------!
188	!
189	numit = numit + nn_fsbc ! Ice model time step
190	!
191	! ! Store previous ice values
192	a_i_b (:,:,:) = a_i (:,:,:) ! ice area
193	e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy
194	v_i_b (:,:,:) = v_i (:,:,:) ! ice volume
195	v_s_b (:,:,:) = v_s (:,:,:) ! snow volume
196	e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy
197	smv_i_b(:,:,:) = smv_i(:,:,:) ! salt content
198	oa_i_b (:,:,:) = oa_i (:,:,:) ! areal age content
199	u_ice_b(:,:) = u_ice(:,:)
200	v_ice_b(:,:) = v_ice(:,:)
201
202	CALL sbc_lim_flx0 ! set diag of mass, heat and salt fluxes to 0
203
204	CALL lim_rst_opn( kt ) ! Open Ice restart file
205	!
206	IF( ln_nicep ) CALL lim_prt( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' ) ! control print
207	! ----------------------------------------------
208	! ice dynamics and transport (except in 1D case)
209	! ----------------------------------------------
210	IF( .NOT. lk_c1d ) THEN
211
212	CALL lim_dyn( kt ) ! Ice dynamics ( rheology/dynamics )
213	CALL lim_trp( kt ) ! Ice transport ( Advection/diffusion )
214	CALL lim_var_glo2eqv ! equivalent variables, requested for rafting
215	IF( ln_nicep ) CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' ) ! control print
216	IF( nn_monocat /= 2 ) &
217	& CALL lim_itd_me ! Mechanical redistribution ! (ridging/rafting)
218	CALL lim_var_agg( 1 )
219	#if defined key_bdy
220	! bdy ice thermo
221	CALL lim_var_glo2eqv ! equivalent variables
222	CALL bdy_ice_lim( kt )
223	CALL lim_var_zapsmall
224	CALL lim_var_agg(1)
225	IF( ln_nicep ) CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' ) ! control print
226	#endif
227	CALL lim_update1
228	ENDIF
229	! !- Change old values for new values
230	u_ice_b(:,:) = u_ice(:,:)
231	v_ice_b(:,:) = v_ice(:,:)
232	a_i_b (:,:,:) = a_i (:,:,:)
233	v_s_b (:,:,:) = v_s (:,:,:)
234	v_i_b (:,:,:) = v_i (:,:,:)
235	e_s_b (:,:,:,:) = e_s (:,:,:,:)
236	e_i_b (:,:,:,:) = e_i (:,:,:,:)
237	oa_i_b (:,:,:) = oa_i (:,:,:)
238	smv_i_b(:,:,:) = smv_i(:,:,:)
239
240	! ----------------------------------------------
241	! ice thermodynamics
242	! ----------------------------------------------
243	CALL lim_var_glo2eqv ! equivalent variables
244	CALL lim_var_agg(1) ! aggregate ice categories
245
246	! previous lead fraction and ice volume for flux calculations
247	pfrld(:,:) = 1._wp - at_i(:,:)
248	phicif(:,:) = vt_i(:,:)
249
250	SELECT CASE( kblk )
251	CASE ( jp_cpl )
252	CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
253	IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice , &
254	& dqns_ice, qla_ice, dqla_ice, nn_limflx )
255	! Latent heat flux is forced to 0 in coupled :
256	! it is included in qns (non-solar heat flux)
257	qla_ice (:,:,:) = 0._wp
258	dqla_ice (:,:,:) = 0._wp
259	END SELECT
260	!
261	CALL lim_var_bv ! bulk brine volume (diag)
262	CALL lim_thd( kt ) ! Ice thermodynamics
263	zcoef = rdt_ice /rday ! Ice natural aging
264	oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
265	IF( ln_nicep ) CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' ) ! control print
266	CALL lim_itd_th( kt ) ! Remap ice categories, lateral accretion
267	CALL lim_update2 ! Global variables update
268
269	IF( ln_nicep ) CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' ) ! control print
270	!
271	CALL lim_sbc_flx( kt ) ! Update surface ocean mass, heat and salt fluxes
272	!
273	IF( ln_nicep ) CALL lim_prt( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' ) ! control print
274	!
275	IF(ln_limdiaout) CALL lim_diahsb ! Diagnostics and outputs
276
277	CALL lim_wri( 1 ) ! Ice outputs
278
279	IF( kt == nit000 .AND. ln_rstart ) &
280	& CALL iom_close( numrir ) ! clem: close input ice restart file
281	!
282	IF( lrst_ice ) CALL lim_rst_write( kt ) ! Ice restart file
283	CALL lim_var_glo2eqv ! ???
284	!
285	IF( ln_nicep ) CALL lim_ctl( kt ) ! alerts in case of model crash
286	!
287	CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
288	!
289	ENDIF ! End sea-ice time step only
290
291	! !--------------------------!
292	! ! at all ocean time step !
293	! !--------------------------!
294	!
295	! ! Update surface ocean stresses (only in ice-dynamic case)
296	! ! otherwise the atm.-ocean stresses are used everywhere
297	IF( ln_limdyn ) CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) ) ! using before instantaneous surf. currents
298	!!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!!
299	!
300	IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim')
301	!
302	END SUBROUTINE sbc_ice_lim
303
304
305	SUBROUTINE sbc_lim_init
306	!!----------------------------------------------------------------------
307	!! * ROUTINE sbc_lim_init *
308	!!
309	!! ** purpose : Allocate all the dynamic arrays of the LIM-3 modules
310	!!----------------------------------------------------------------------
311	INTEGER :: ierr
312	!!----------------------------------------------------------------------
313	IF(lwp) WRITE(numout,*)
314	IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'
315	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ via Louvain la Neuve Ice Model (LIM-3) time stepping'
316	!
317	!
318	! ! Allocate the ice arrays
319	ierr = ice_alloc () ! ice variables
320	ierr = ierr + dom_ice_alloc () ! domain
321	ierr = ierr + sbc_ice_alloc () ! surface forcing
322	ierr = ierr + thd_ice_alloc () ! thermodynamics
323	ierr = ierr + lim_itd_me_alloc () ! ice thickness distribution - mechanics
324	!
325	IF( lk_mpp ) CALL mpp_sum( ierr )
326	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'sbc_lim_init : unable to allocate ice arrays')
327	!
328	! ! adequation jpk versus ice/snow layers/categories
329	IF( jpl > jpk .OR. (nlay_i+1) > jpk .OR. nlay_s > jpk ) &
330	& CALL ctl_stop( 'STOP', &
331	& 'sbc_lim_init: the 3rd dimension of workspace arrays is too small.', &
332	& 'use more ocean levels or less ice/snow layers/categories.' )
333
334	! Open the reference and configuration namelist files and namelist output file
335	CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
336	CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
337	IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
338	!
339	CALL ice_run ! set some ice run parameters
340	!
341	CALL lim_thd_init ! set ice thermodynics parameters
342	!
343	CALL lim_thd_sal_init ! set ice salinity parameters
344	!
345	rdt_ice = nn_fsbc * rdttra(1) ! sea-ice timestep
346	r1_rdtice = 1._wp / rdt_ice ! sea-ice timestep inverse
347	!
348	CALL lim_msh ! ice mesh initialization
349	!
350	CALL lim_itd_init ! ice thickness distribution initialization
351	!
352	CALL lim_itd_me_init ! ice thickness distribution initialization
353	! ! Initial sea-ice state
354	IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst
355	numit = 0
356	numit = nit000 - 1
357	CALL lim_istate
358	CALL lim_var_agg(1)
359	CALL lim_var_glo2eqv
360	ELSE ! start from a restart file
361	CALL lim_rst_read
362	numit = nit000 - 1
363	CALL lim_var_agg(1)
364	CALL lim_var_glo2eqv
365	ENDIF
366	!
367	CALL lim_sbc_init ! ice surface boundary condition
368	!
369	fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction
370	tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu
371	!
372	nstart = numit + nn_fsbc
373	nitrun = nitend - nit000 + 1
374	nlast = numit + nitrun
375	!
376	IF( nstock == 0 ) nstock = nlast + 1
377	!
378	END SUBROUTINE sbc_lim_init
379
380
381	SUBROUTINE ice_run
382	!!-------------------------------------------------------------------
383	!! * ROUTINE ice_run *
384	!!
385	!! ** Purpose : Definition some run parameter for ice model
386	!!
387	!! ** Method : Read the namicerun namelist and check the parameter
388	!! values called at the first timestep (nit000)
389	!!
390	!! ** input : Namelist namicerun
391	!!-------------------------------------------------------------------
392	NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout
393	INTEGER :: ios ! Local integer output status for namelist read
394	!!-------------------------------------------------------------------
395	!
396	REWIND( numnam_ice_ref ) ! Namelist namicerun in reference namelist : Parameters for ice
397	READ ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
398	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
399
400	REWIND( numnam_ice_cfg ) ! Namelist namicerun in configuration namelist : Parameters for ice
401	READ ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
402	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
403	IF(lwm) WRITE ( numoni, namicerun )
404	!
405	!
406	IF(lwp) THEN ! control print
407	WRITE(numout,*)
408	WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
409	WRITE(numout,*) ' ~~~~~~'
410	WRITE(numout,*) ' switch for ice dynamics (1) or not (0) ln_limdyn = ', ln_limdyn
411	WRITE(numout,*) ' maximum ice concentration = ', amax
412	WRITE(numout,*) ' Several ice points in the ice or not in ocean.output = ', ln_nicep
413	WRITE(numout,*) ' Diagnose heat/salt budget or not ln_limdiahsb = ', ln_limdiahsb
414	WRITE(numout,*) ' Output heat/salt budget or not ln_limdiaout = ', ln_limdiaout
415	ENDIF
416	!
417	!IF( lk_mpp .AND. ln_nicep ) THEN
418	! ln_nicep = .FALSE.
419	! CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' )
420	!ENDIF
421	IF( ln_nicep ) THEN ! control print at a given point
422	jiindx = 15 ; jjindx = 44
423	IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
424	ENDIF
425	!
426	END SUBROUTINE ice_run
427
428
429	SUBROUTINE lim_itd_init
430	!!------------------------------------------------------------------
431	!! * ROUTINE lim_itd_init *
432	!!
433	!! ** Purpose : Initializes the ice thickness distribution
434	!! ** Method : ...
435	!! Note : hi_max(jpl) is here set up to a value close to 7 m for
436	!! limistate (only) and is changed to 99 m in sbc_lim_init
437	!!------------------------------------------------------------------
438	INTEGER :: jl ! dummy loop index
439	REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars
440	REAL(wp) :: zhmax, znum, zden, zalpha !
441	!!------------------------------------------------------------------
442
443	IF(lwp) WRITE(numout,*)
444	IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice thickness distribution '
445	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
446
447	!------------------------------------------------------------------------------!
448	! 1) Ice thickness distribution parameters initialization
449	!------------------------------------------------------------------------------!
450	IF(lwp) WRITE(numout,*) ' Number of ice categories jpl = ', jpl
451
452	!- Thickness categories boundaries
453	!----------------------------------
454	! Clem - je sais pas encore dans quelle namelist les mettre, ca depend des chgts liés à bdy
455	nn_hibnd = 2 ! thickness category boundaries: tanh function (1) h^(-alpha) (2)
456	rn_hibnd = 2.5 ! mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only)
457
458	hi_max(:) = 0._wp
459
460	SELECT CASE ( nn_hibnd )
461	!----------------------
462	CASE (1) ! tanh function (CICE)
463	!----------------------
464	zc1 = 3._wp / REAL( jpl, wp )
465	zc2 = 10._wp * zc1
466	zc3 = 3._wp
467
468	DO jl = 1, jpl
469	zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
470	hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
471	END DO
472
473	!----------------------
474	CASE (2) ! h^(-alpha) function
475	!----------------------
476	zalpha = 0.05 ! exponent of the transform function
477
478	zhmax = 3.*rn_hibnd
479
480	DO jl = 1, jpl
481	znum = jpl * ( zhmax+1 )**zalpha
482	zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl
483	hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
484	END DO
485
486	END SELECT
487
488	DO jl = 1, jpl
489	hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
490	END DO
491	! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl)
492	hi_max(jpl) = 99._wp
493
494	IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
495	IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
496	!
497	END SUBROUTINE lim_itd_init
498
499
500	SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, &
501	& pdqn_ice, pqla_ice, pdql_ice, k_limflx )
502	!!---------------------------------------------------------------------
503	!! * ROUTINE sbc_ice_lim *
504	!!
505	!! ** Purpose : update the ice surface boundary condition by averaging and / or
506	!! redistributing fluxes on ice categories
507	!!
508	!! ** Method : average then redistribute
509	!!
510	!! ** Action :
511	!!---------------------------------------------------------------------
512	INTEGER , INTENT(in ) :: k_limflx ! =-1 do nothing; =0 average ;
513	! =1 average and redistribute ; =2 redistribute
514	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature
515	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo
516	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqns_ice ! non solar flux
517	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux
518	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity
519	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqla_ice ! latent heat flux
520	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdql_ice ! latent heat flux sensitivity
521	!
522	INTEGER :: jl ! dummy loop index
523	!
524	REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m ! Mean albedo over all categories
525	REAL(wp), POINTER, DIMENSION(:,:) :: ztem_m ! Mean temperature over all categories
526	!
527	REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories
528	REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m ! Mean non solar heat flux over all categories
529	REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_m ! Mean latent heat flux over all categories
530	REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m ! Mean d(qns)/dT over all categories
531	REAL(wp), POINTER, DIMENSION(:,:) :: z_dql_m ! Mean d(qla)/dT over all categories
532	!!----------------------------------------------------------------------
533
534	IF( nn_timing == 1 ) CALL timing_start('ice_lim_flx')
535	!
536	!
537	SELECT CASE( k_limflx ) !== averaged on all ice categories ==!
538	CASE( 0 , 1 )
539	CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m)
540	!
541	z_qns_m(:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
542	z_qsr_m(:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
543	z_dqn_m(:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
544	z_qla_m(:,:) = fice_ice_ave ( pqla_ice (:,:,:) )
545	z_dql_m(:,:) = fice_ice_ave ( pdql_ice (:,:,:) )
546	DO jl = 1, jpl
547	pdqn_ice(:,:,jl) = z_dqn_m(:,:)
548	pdql_ice(:,:,jl) = z_dql_m(:,:)
549	END DO
550	!
551	DO jl = 1, jpl
552	pqns_ice(:,:,jl) = z_qns_m(:,:)
553	pqsr_ice(:,:,jl) = z_qsr_m(:,:)
554	pqla_ice(:,:,jl) = z_qla_m(:,:)
555	END DO
556	!
557	CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m)
558	END SELECT
559
560	SELECT CASE( k_limflx ) !== redistribution on all ice categories ==!
561	CASE( 1 , 2 )
562	CALL wrk_alloc( jpi,jpj, zalb_m, ztem_m )
563	!
564	zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )
565	ztem_m(:,:) = fice_ice_ave ( ptn_ice (:,:,:) )
566	DO jl = 1, jpl
567	pqns_ice(:,:,jl) = pqns_ice(:,:,jl) + pdqn_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
568	pqla_ice(:,:,jl) = pqla_ice(:,:,jl) + pdql_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
569	pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
570	END DO
571	!
572	CALL wrk_dealloc( jpi,jpj, zalb_m, ztem_m )
573	END SELECT
574	!
575	IF( nn_timing == 1 ) CALL timing_stop('ice_lim_flx')
576	!
577	END SUBROUTINE ice_lim_flx
578
579	SUBROUTINE sbc_lim_flx0
580	!!----------------------------------------------------------------------
581	!! * ROUTINE sbc_lim_flx0 *
582	!!
583	!! ** purpose : set ice-ocean and ice-atm. fluxes to zeros at the beggining
584	!! of the time step
585	!!----------------------------------------------------------------------
586	sfx (:,:) = 0._wp ;
587	sfx_bri(:,:) = 0._wp ;
588	sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp
589	sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp
590	sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp
591	sfx_res(:,:) = 0._wp
592
593	wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp
594	wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp
595	wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp
596	wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp
597	wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp
598	wfx_spr(:,:) = 0._wp ;
599
600	hfx_in (:,:) = 0._wp ; hfx_out(:,:) = 0._wp
601	hfx_thd(:,:) = 0._wp ;
602	hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp
603	hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp
604	hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp
605	hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp
606	hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp
607	hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp
608
609	afx_tot(:,:) = 0._wp ;
610	afx_dyn(:,:) = 0._wp ; afx_thd(:,:) = 0._wp
611
612	END SUBROUTINE sbc_lim_flx0
613
614	FUNCTION fice_cell_ave ( ptab )
615	!!--------------------------------------------------------------------------
616	!! * Compute average over categories, for grid cell (ice covered and free ocean)
617	!!--------------------------------------------------------------------------
618	REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
619	REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
620	INTEGER :: jl ! Dummy loop index
621
622	fice_cell_ave (:,:) = 0.0_wp
623
624	DO jl = 1, jpl
625	fice_cell_ave (:,:) = fice_cell_ave (:,:) &
626	& + a_i (:,:,jl) * ptab (:,:,jl)
627	END DO
628
629	END FUNCTION fice_cell_ave
630
631
632	FUNCTION fice_ice_ave ( ptab )
633	!!--------------------------------------------------------------------------
634	!! * Compute average over categories, for ice covered part of grid cell
635	!!--------------------------------------------------------------------------
636	REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
637	REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
638
639	fice_ice_ave (:,:) = 0.0_wp
640	WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
641
642	END FUNCTION fice_ice_ave
643
644
645	#else
646	!!----------------------------------------------------------------------
647	!! Default option Dummy module NO LIM 3.0 sea-ice model
648	!!----------------------------------------------------------------------
649	CONTAINS
650	SUBROUTINE sbc_ice_lim ( kt, kblk ) ! Dummy routine
651	WRITE(,) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk
652	END SUBROUTINE sbc_ice_lim
653	#endif
654
655	!!======================================================================
656	END MODULE sbcice_lim

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

Download in other formats: