Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbcice_lim.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 33.5 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 7910

Download in other formats: