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.
icestp.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90 @ 8410

Last change on this file since 8410 was 8410, checked in by clem, 7 years ago

correction on the last commit
File size: 28.1 KB
Line 
1MODULE icestp
2   !!======================================================================
3   !!                       ***  MODULE  icestp  ***
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   !!   ice_stp  : 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 iceforcing      ! Surface boundary condition for sea ice
31   !
32   USE phycst          ! Define parameters for the routines
33   USE eosbn2          ! equation of state
34   USE icerhg          ! Ice rheology
35   USE iceadv          ! Ice advection
36   USE limthd          ! Ice thermodynamics
37   USE icerdgrft       ! Ice ridging/rafting
38   USE limsbc          ! sea surface boundary condition
39   USE limdiahsb       ! Ice budget diagnostics
40   USE limwri          ! Ice outputs
41   USE limrst          ! Ice restarts
42   USE iceerr1         ! Ice corrections after dynamics
43   USE iceerr2         ! Ice corrections after thermo
44   USE limvar          ! Ice variables switch
45   USE limctl          !
46   ! MV MP 2016
47   USE limmp
48   ! END MV MP 2016
49   USE limistate       ! LIM initial state
50   USE limthd_sal      ! LIM ice thermodynamics: salinity
51   !
52   USE c1d             ! 1D vertical configuration
53   USE in_out_manager  ! I/O manager
54   USE iom             ! I/O manager library
55   USE prtctl          ! Print control
56   USE lib_fortran     !
57   USE lbclnk          ! lateral boundary condition - MPP link
58   USE lib_mpp         ! MPP library
59   USE timing          ! Timing
60
61   USE bdy_oce   , ONLY: ln_bdy
62   USE bdyice          ! unstructured open boundary data
63# if defined key_agrif
64   USE agrif_ice
65   USE agrif_lim3_update
66   USE agrif_lim3_interp
67# endif
68
69   IMPLICIT NONE
70   PRIVATE
71
72   PUBLIC ice_stp  ! routine called by sbcmod.F90
73   PUBLIC ice_init ! routine called by sbcmod.F90
74
75   !! * Substitutions
76#  include "vectopt_loop_substitute.h90"
77   !!----------------------------------------------------------------------
78   !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
79   !! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $
80   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
81   !!----------------------------------------------------------------------
82CONTAINS
83
84   SUBROUTINE ice_stp( kt, ksbc )
85      !!---------------------------------------------------------------------
86      !!                  ***  ROUTINE ice_stp  ***
87      !!
88      !! ** Purpose :   update the ocean surface boundary condition via the
89      !!                Louvain la Neuve Sea Ice Model time stepping
90      !!
91      !! ** Method  :   ice model time stepping
92      !!              - call the ice dynamics routine
93      !!              - call the ice advection/diffusion routine
94      !!              - call the ice thermodynamics routine
95      !!              - call the routine that computes mass and
96      !!                heat fluxes at the ice/ocean interface
97      !!              - save the outputs
98      !!              - save the outputs for restart when necessary
99      !!
100      !! ** Action  : - time evolution of the LIM sea-ice model
101      !!              - update all sbc variables below sea-ice:
102      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx
103      !!---------------------------------------------------------------------
104      INTEGER, INTENT(in) ::   kt      ! ocean time step
105      INTEGER, INTENT(in) ::   ksbc    ! type of sbc flux ( 1 = user defined formulation,
106                                       !                    3 = bulk formulation,
107                                       !                    4 = Pure Coupled formulation)
108      INTEGER  ::   jl                 ! dummy loop index
109      !!----------------------------------------------------------------------
110
111      IF( nn_timing == 1 )   CALL timing_start('ice_stp')
112
113      !-----------------------!
114      ! --- Ice time step --- !
115      !-----------------------!
116      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
117         
118         ! Ice model time step
119         kt_ice = kt
120
121# if defined key_agrif
122         IF( .NOT. Agrif_Root() )  lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
123# endif
124
125         ! mean surface ocean current at ice velocity point (C-grid dynamics :  U- & V-points as the ocean)
126         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
127         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
128
129         ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
130         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
131         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
132
133         ! Mask sea ice surface temperature (set to rt0 over land)
134         DO jl = 1, jpl
135            t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
136         END DO
137         !
138                                      CALL ice_bef         ! Store previous ice values
139         !------------------------------------------------!
140         ! --- Dynamical coupling with the atmosphere --- !
141         !------------------------------------------------!
142         ! it provides:
143         !    utau_ice, vtau_ice = surface ice stress [N/m2]
144         !--------------------------------------------------
145                                      CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )
146                                     
147         !-------------------------------------------------------!
148         ! --- ice dynamics and transport (except in 1D case) ---!
149         !-------------------------------------------------------!
150                                      CALL ice_diag0           ! set diag of mass, heat and salt fluxes to 0
151                                      CALL lim_rst_opn( kt )   ! Open Ice restart file
152         !
153         ! --- zap this if no ice dynamics --- !
154         IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN
155                                      CALL ice_rhg( kt )       ! -- rheology 
156                                      CALL ice_adv( kt )       ! -- advection
157            IF( nn_limdyn == 2 .AND. nn_monocat /= 2 )  &      ! -- ridging/rafting
158               &                      CALL ice_rdgrft         
159            IF( nn_limdyn == 2 )      CALL ice_err1( kt )      ! -- Corrections
160            !
161         ENDIF
162
163         ! ---
164#if defined key_agrif
165         IF( .NOT. Agrif_Root() )     CALL agrif_interp_lim3('T')
166#endif
167         IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice( kt )       ! -- bdy ice thermo
168         ! previous lead fraction and ice volume for flux calculations
169                                      CALL lim_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation
170                                      CALL lim_var_agg(1)      ! at_i for coupling
171                                      CALL ice_bef
172
173         !------------------------------------------------------!
174         ! --- Thermodynamical coupling with the atmosphere --- !
175         !------------------------------------------------------!
176         ! It provides the following fields used in sea ice model:
177         !    fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%]
178         !    emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
179         !    sprecip                                  = solid precipitation                           [Kg/m2/s]
180         !    evap_ice                                 = sublimation                                   [Kg/m2/s]
181         !    qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2]
182         !    qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2]
183         !    dqns_ice                                 = non solar  heat sensistivity                  [W/m2]
184         !    qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
185         !------------------------------------------------------------------------------------------------------
186                                      CALL ice_forcing_flx( kt, ksbc )
187
188         !----------------------------!
189         ! --- ice thermodynamics --- !
190         !----------------------------!
191         ! --- zap this if no ice thermo --- !
192         IF( ln_limthd )              CALL lim_thd( kt )        ! -- Ice thermodynamics     
193
194         ! MV MP 2016
195         IF ( ln_pnd )                CALL lim_mp( kt )         ! -- Melt ponds
196         ! END MV MP 2016
197
198         IF( ln_limthd )              CALL ice_err2( kt )       ! -- Corrections
199         ! ---
200# if defined key_agrif
201         IF( .NOT. Agrif_Root() )     CALL agrif_update_lim3( kt )
202# endif
203                                      CALL lim_var_glo2eqv      ! necessary calls (at least for coupling)
204                                      CALL lim_var_agg( 2 )     ! necessary calls (at least for coupling)
205                                      !
206!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
207!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
208!!# if defined key_agrif
209!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid()
210!!# endif
211                                      CALL lim_sbc_flx( kt )    ! -- Update surface ocean mass, heat and salt fluxes
212!!# if defined key_agrif
213!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
214!!# endif
215         IF( ln_limdiahsb )           CALL lim_diahsb( kt )     ! -- Diagnostics and outputs
216         !
217                                      CALL lim_wri( 1 )         ! -- Ice outputs
218         !
219         IF( kt == nit000 .AND. ln_rstart )   &
220            &                         CALL iom_close( numrir )  ! close input ice restart file
221         !
222         IF( lrst_ice )               CALL lim_rst_write( kt )  ! -- Ice restart file
223         !
224         IF( ln_limctl )              CALL lim_ctl( kt )        ! alerts in case of model crash
225         !
226      ENDIF   ! End sea-ice time step only
227
228      !-------------------------!
229      ! --- Ocean time step --- !
230      !-------------------------!
231      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
232      !    using before instantaneous surf. currents
233      IF( ln_limdyn )                 CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )
234!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
235      !
236      IF( nn_timing == 1 ) CALL timing_stop('ice_stp')
237      !
238   END SUBROUTINE ice_stp
239
240
241   SUBROUTINE ice_init
242      !!----------------------------------------------------------------------
243      !!                  ***  ROUTINE ice_init  ***
244      !!
245      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules
246      !!----------------------------------------------------------------------
247      INTEGER :: ji, jj, ierr
248      !!----------------------------------------------------------------------
249      IF(lwp) WRITE(numout,*)
250      IF(lwp) WRITE(numout,*) 'ice_init : update ocean surface boudary condition' 
251      IF(lwp) WRITE(numout,*) '~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
252      !
253      !                                ! Open the reference and configuration namelist files and namelist output file
254      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
255      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
256      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
257      !
258      CALL ice_run_init                ! set some ice run parameters
259      !
260      !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
261      ierr =        ice_alloc        ()      ! ice variables
262      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
263      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics
264      IF( ln_limdyn )   ierr = ierr + ice_rdgrft_alloc ()      ! ridging/rafting
265      !
266      IF( lk_mpp    )   CALL mpp_sum( ierr )
267      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
268      !
269      CALL ice_rhg_init                ! set ice dynamics parameters
270      !
271      CALL ice_itd_init                ! ice thickness distribution initialization
272      !
273      CALL lim_thd_init                ! set ice thermodynics parameters
274      !
275      CALL lim_thd_sal_init            ! set ice salinity parameters
276       
277      ! MV MP 2016
278      CALL lim_mp_init                 ! set melt ponds parameters
279      ! END MV MP 2016
280
281      IF( ln_limdyn )   CALL ice_rdgrft_init             ! ice thickness distribution initialization for ridging/rafting
282      !                                ! Initial sea-ice state
283      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
284         CALL lim_istate
285      ELSE                                    ! start from a restart file
286         CALL lim_rst_read
287      ENDIF
288      CALL lim_var_agg(2)
289      CALL lim_var_glo2eqv
290      !
291      CALL lim_sbc_init                 ! ice surface boundary condition
292      !
293      IF( ln_limdiahsb) CALL lim_diahsb_init  ! initialization for diags
294      !
295      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction
296      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
297      !
298      DO jj = 1, jpj
299         DO ji = 1, jpi
300            IF( gphit(ji,jj) > 0._wp ) THEN  ;  rn_amax_2d(ji,jj) = rn_amax_n  ! NH
301            ELSE                             ;  rn_amax_2d(ji,jj) = rn_amax_s  ! SH
302            ENDIF
303         END DO
304      END DO
305      !
306   END SUBROUTINE ice_init
307
308
309   SUBROUTINE ice_run_init
310      !!-------------------------------------------------------------------
311      !!                  ***  ROUTINE ice_run_init ***
312      !!
313      !! ** Purpose :   Definition some run parameter for ice model
314      !!
315      !! ** Method  :   Read the namicerun namelist and check the parameter
316      !!              values called at the first timestep (nit000)
317      !!
318      !! ** input   :   Namelist namicerun
319      !!-------------------------------------------------------------------
320      INTEGER  ::   ios                 ! Local integer output status for namelist read
321      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   &
322         &                cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice 
323      NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt 
324      !!-------------------------------------------------------------------
325      !
326      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice
327      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
328901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
329
330      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice
331      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
332902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
333      IF(lwm) WRITE ( numoni, namicerun )
334      !
335      REWIND( numnam_ice_ref )              ! Namelist namicediag in reference namelist : Parameters for ice
336      READ  ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
337903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
338
339      REWIND( numnam_ice_cfg )              ! Namelist namicediag in configuration namelist : Parameters for ice
340      READ  ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
341904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
342      IF(lwm) WRITE ( numoni, namicediag )
343      !
344      IF(lwp) THEN                        ! control print
345         WRITE(numout,*)
346         WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
347         WRITE(numout,*) ' ~~~~~~'
348         WRITE(numout,*) '   number of ice  categories                              jpl    = ', jpl
349         WRITE(numout,*) '   number of ice  layers                                  nlay_i = ', nlay_i
350         WRITE(numout,*) '   number of snow layers                                  nlay_s = ', nlay_s
351         WRITE(numout,*) '   virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat
352         WRITE(numout,*) '   maximum ice concentration for NH                              = ', rn_amax_n 
353         WRITE(numout,*) '   maximum ice concentration for SH                              = ', rn_amax_s
354         WRITE(numout,*) '   Ice thermodynamics (T) or not (F)                  ln_limthd  = ', ln_limthd
355         WRITE(numout,*) '   Ice dynamics       (T) or not (F)                  ln_limdyn  = ', ln_limdyn
356         WRITE(numout,*) '     (ln_limdyn=T) Ice dynamics switch                nn_limdyn  = ', nn_limdyn
357         WRITE(numout,*) '       2: total'
358         WRITE(numout,*) '       1: advection only (no diffusion, no ridging/rafting)'
359         WRITE(numout,*) '       0: advection only (as 1 + prescribed velocity, bypass rheology)'
360         WRITE(numout,*) '     (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0)           = ', rn_uice
361         WRITE(numout,*) '     (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0)           = ', rn_vice
362         WRITE(numout,*)
363         WRITE(numout,*) '...and ice diagnostics'
364         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
365         WRITE(numout,*) '   Diagnose online heat/mass/salt budget     ln_limdiachk  = ', ln_limdiachk
366         WRITE(numout,*) '   Output          heat/mass/salt budget     ln_limdiahsb  = ', ln_limdiahsb
367         WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl
368         WRITE(numout,*) '   i-index for control prints (ln_limctl=true)             = ', iiceprt
369         WRITE(numout,*) '   j-index for control prints (ln_limctl=true)             = ', jiceprt
370      ENDIF
371      !
372      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
373         nn_monocat = 0
374         IF(lwp) WRITE(numout,*)
375         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
376      ENDIF
377      IF ( ( jpl == 1 ) .AND. ( nn_monocat == 0 ) ) THEN
378         CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
379      ENDIF
380      !
381      ! sea-ice timestep and inverse
382      rdt_ice   = REAL(nn_fsbc) * rdt 
383      r1_rdtice = 1._wp / rdt_ice 
384
385      ! inverse of nlay_i and nlay_s
386      r1_nlay_i = 1._wp / REAL( nlay_i, wp )
387      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
388      !
389      IF( lwp .AND. ln_bdy .AND. ln_limdiachk )  &
390         &   CALL ctl_warn('online conservation check activated but it does not work with BDY')
391      !
392      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice
393      !
394   END SUBROUTINE ice_run_init
395
396
397   SUBROUTINE ice_itd_init
398      !!------------------------------------------------------------------
399      !!                ***  ROUTINE ice_itd_init ***
400      !!
401      !! ** Purpose :   Initializes the ice thickness distribution
402      !! ** Method  :   ...
403      !! ** input   :   Namelist namiceitd
404      !!-------------------------------------------------------------------
405      INTEGER  ::   ios                 ! Local integer output status for namelist read
406      NAMELIST/namiceitd/ rn_himean
407      !
408      INTEGER  ::   jl                   ! dummy loop index
409      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
410      REAL(wp) ::   zhmax, znum, zden, zalpha !
411      !!------------------------------------------------------------------
412      !
413      REWIND( numnam_ice_ref )              ! Namelist namiceitd in reference namelist : Parameters for ice
414      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
415905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
416
417      REWIND( numnam_ice_cfg )              ! Namelist namiceitd in configuration namelist : Parameters for ice
418      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
419906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
420      IF(lwm) WRITE ( numoni, namiceitd )
421      !
422      IF(lwp) THEN                        ! control print
423         WRITE(numout,*)
424         WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution '
425         WRITE(numout,*) '~~~~~~~~~~~~'
426         WRITE(numout,*) '   mean ice thickness in the domain            rn_himean = ', rn_himean
427      ENDIF
428      !
429      !----------------------------------
430      !- Thickness categories boundaries
431      !----------------------------------
432      !
433      !==  h^(-alpha) function  ==!
434      zalpha = 0.05_wp
435      zhmax  = 3._wp * rn_himean
436      DO jl = 1, jpl
437         znum = jpl * ( zhmax+1 )**zalpha
438         zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
439         hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
440      END DO
441      !
442      DO jl = 1, jpl                ! mean thickness by category
443         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
444      END DO
445      !
446      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
447      !
448      IF(lwp) WRITE(numout,*)
449      IF(lwp) WRITE(numout,*) '      Thickness category boundaries '
450      IF(lwp) WRITE(numout,*) '         hi_max ', hi_max(0:jpl)
451      !
452   END SUBROUTINE ice_itd_init
453
454   SUBROUTINE ice_bef
455      !!----------------------------------------------------------------------
456      !!                  ***  ROUTINE ice_bef  ***
457      !!
458      !! ** purpose :  store ice variables at "before" time step
459      !!----------------------------------------------------------------------
460      INTEGER  ::   ji, jj, jl      ! dummy loop index
461      !!----------------------------------------------------------------------
462      !
463      DO jl = 1, jpl
464         DO jj = 2, jpjm1
465            DO ji = 2, jpim1
466               a_i_b  (ji,jj,jl)   = a_i  (ji,jj,jl)     ! ice area
467               v_i_b  (ji,jj,jl)   = v_i  (ji,jj,jl)     ! ice volume
468               v_s_b  (ji,jj,jl)   = v_s  (ji,jj,jl)     ! snow volume
469               smv_i_b(ji,jj,jl)   = smv_i(ji,jj,jl)     ! salt content
470               oa_i_b (ji,jj,jl)   = oa_i (ji,jj,jl)     ! areal age content
471               e_s_b  (ji,jj,:,jl) = e_s  (ji,jj,:,jl)   ! snow thermal energy
472               e_i_b  (ji,jj,:,jl) = e_i  (ji,jj,:,jl)   ! ice thermal energy
473               !                                         ! ice thickness
474               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes
475               ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
476               ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
477            END DO
478         END DO   
479      END DO
480      CALL lbc_lnk_multi( a_i_b, 'T', 1., v_i_b , 'T', 1., v_s_b , 'T', 1., smv_i_b, 'T', 1., &
481         &               oa_i_b, 'T', 1., ht_i_b, 'T', 1., ht_s_b, 'T', 1. )
482      CALL lbc_lnk( e_i_b, 'T', 1. )
483      CALL lbc_lnk( e_s_b, 'T', 1. )
484     
485      ! ice velocities & total concentration
486      DO jj = 2, jpjm1
487         DO ji = 2, jpim1
488            at_i_b(ji,jj)  = SUM( a_i_b(ji,jj,:) )
489            u_ice_b(ji,jj) = u_ice(ji,jj)
490            v_ice_b(ji,jj) = v_ice(ji,jj)
491         END DO
492      END DO
493      CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. )
494     
495   END SUBROUTINE ice_bef
496
497
498   SUBROUTINE ice_diag0
499      !!----------------------------------------------------------------------
500      !!                  ***  ROUTINE ice_diag0  ***
501      !!
502      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
503      !!               of the time step
504      !!----------------------------------------------------------------------
505      INTEGER  ::   ji, jj      ! dummy loop index
506      !!----------------------------------------------------------------------
507      DO jj = 1, jpj
508         DO ji = 1, jpi
509            sfx    (ji,jj) = 0._wp   ;
510            sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp
511            sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp
512            sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp
513            sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp
514            sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp
515            !
516            wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp
517            wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp
518            wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp
519            wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp
520            wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp
521            wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp 
522            wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
523            wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
524            wfx_snw_sni(ji,jj) = 0._wp 
525            ! MV MP 2016
526            wfx_pnd(ji,jj) = 0._wp
527            ! END MV MP 2016
528           
529            hfx_thd(ji,jj) = 0._wp   ;
530            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
531            hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp
532            hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp
533            hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp
534            hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp
535            hfx_err(ji,jj) = 0._wp   ;   hfx_err_rem(ji,jj) = 0._wp
536            hfx_err_dif(ji,jj) = 0._wp
537            wfx_err_sub(ji,jj) = 0._wp
538            !
539            afx_tot(ji,jj) = 0._wp   ;
540            afx_dyn(ji,jj) = 0._wp   ;   afx_thd(ji,jj) = 0._wp
541            !
542            diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp
543            diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp
544           
545            ! SIMIP diagnostics
546            diag_dms_dyn(ji,jj)  = 0._wp ; diag_dmi_dyn(ji,jj)  = 0._wp
547            diag_fc_bo(ji,jj)    = 0._wp ; diag_fc_su(ji,jj)    = 0._wp
548           
549            tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
550         END DO
551      END DO
552     
553   END SUBROUTINE ice_diag0
554
555#else
556   !!----------------------------------------------------------------------
557   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
558   !!----------------------------------------------------------------------
559CONTAINS
560   !
561   SUBROUTINE ice_stp ( kt, ksbc )     ! Dummy routine
562      INTEGER, INTENT(in) ::   kt, ksbc
563      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
564   END SUBROUTINE ice_stp
565   !
566   SUBROUTINE ice_init                 ! Dummy routine
567   END SUBROUTINE ice_init
568#endif
569
570   !!======================================================================
571END MODULE icestp
Note: See TracBrowser for help on using the repository browser.