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

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

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

Download in other formats: