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

Last change on this file since 8426 was 8426, checked in by clem, 3 years ago

last routine names to be changed

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