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

Last change on this file since 8407 was 8407, checked in by clem, 5 years ago

change calls in icestp.F90 for rheology

File size: 28.6 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 dynamics
35   USE limtrp          ! Ice transport
36   USE limthd          ! Ice thermodynamics
37   USE limitd_me       ! Mechanics on ice thickness distribution
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 limupdate1      ! update of global variables
43   USE limupdate2      ! update of global variables
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            !
156            IF( nn_limdyn /= 0 ) THEN                          ! -- Ice dynamics
157                                      CALL ice_rhg( kt )       !     rheology 
158            ELSE
159               u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity
160               v_ice(:,:) = rn_vice * vmask(:,:,1)
161               !!CALL RANDOM_NUMBER(u_ice(:,:))
162               !!CALL RANDOM_NUMBER(v_ice(:,:))
163            ENDIF
164                                      CALL lim_trp( kt )       ! -- Ice transport (Advection/diffusion)
165            IF( nn_limdyn == 2 .AND. nn_monocat /= 2 )  &      ! -- Mechanical redistribution (ridging/rafting)
166               &                      CALL lim_itd_me         
167            IF( nn_limdyn == 2 )      CALL lim_update1( kt )   ! -- Corrections
168            !
169         ENDIF
170
171         ! ---
172#if defined key_agrif
173         IF( .NOT. Agrif_Root() )     CALL agrif_interp_lim3('T')
174#endif
175         IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice( kt )       ! -- bdy ice thermo
176         ! previous lead fraction and ice volume for flux calculations
177                                      CALL lim_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation
178                                      CALL lim_var_agg(1)      ! at_i for coupling
179                                      CALL ice_bef
180
181         !------------------------------------------------------!
182         ! --- Thermodynamical coupling with the atmosphere --- !
183         !------------------------------------------------------!
184         ! It provides the following fields used in sea ice model:
185         !    fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%]
186         !    emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s]
187         !    sprecip                                  = solid precipitation                           [Kg/m2/s]
188         !    evap_ice                                 = sublimation                                   [Kg/m2/s]
189         !    qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2]
190         !    qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2]
191         !    dqns_ice                                 = non solar  heat sensistivity                  [W/m2]
192         !    qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
193         !------------------------------------------------------------------------------------------------------
194                                      CALL ice_forcing_flx( kt, ksbc )
195
196         !----------------------------!
197         ! --- ice thermodynamics --- !
198         !----------------------------!
199         ! --- zap this if no ice thermo --- !
200         IF( ln_limthd )              CALL lim_thd( kt )        ! -- Ice thermodynamics     
201
202         ! MV MP 2016
203         IF ( ln_pnd )                CALL lim_mp( kt )         ! -- Melt ponds
204         ! END MV MP 2016
205
206         IF( ln_limthd )              CALL lim_update2( kt )    ! -- Corrections
207         ! ---
208# if defined key_agrif
209         IF( .NOT. Agrif_Root() )     CALL agrif_update_lim3( kt )
210# endif
211                                      CALL lim_var_glo2eqv      ! necessary calls (at least for coupling)
212                                      CALL lim_var_agg( 2 )     ! necessary calls (at least for coupling)
213                                      !
214!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
215!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
216!!# if defined key_agrif
217!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid()
218!!# endif
219                                      CALL lim_sbc_flx( kt )    ! -- Update surface ocean mass, heat and salt fluxes
220!!# if defined key_agrif
221!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
222!!# endif
223         IF( ln_limdiahsb )           CALL lim_diahsb( kt )     ! -- Diagnostics and outputs
224         !
225                                      CALL lim_wri( 1 )         ! -- Ice outputs
226         !
227         IF( kt == nit000 .AND. ln_rstart )   &
228            &                         CALL iom_close( numrir )  ! close input ice restart file
229         !
230         IF( lrst_ice )               CALL lim_rst_write( kt )  ! -- Ice restart file
231         !
232         IF( ln_limctl )              CALL lim_ctl( kt )        ! alerts in case of model crash
233         !
234      ENDIF   ! End sea-ice time step only
235
236      !-------------------------!
237      ! --- Ocean time step --- !
238      !-------------------------!
239      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
240      !    using before instantaneous surf. currents
241      IF( ln_limdyn )                 CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )
242!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
243      !
244      IF( nn_timing == 1 ) CALL timing_stop('ice_stp')
245      !
246   END SUBROUTINE ice_stp
247
248
249   SUBROUTINE ice_init
250      !!----------------------------------------------------------------------
251      !!                  ***  ROUTINE ice_init  ***
252      !!
253      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules
254      !!----------------------------------------------------------------------
255      INTEGER :: ji, jj, ierr
256      !!----------------------------------------------------------------------
257      IF(lwp) WRITE(numout,*)
258      IF(lwp) WRITE(numout,*) 'ice_init : update ocean surface boudary condition' 
259      IF(lwp) WRITE(numout,*) '~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
260      !
261      !                                ! Open the reference and configuration namelist files and namelist output file
262      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
263      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
264      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
265      !
266      CALL ice_run_init                ! set some ice run parameters
267      !
268      !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
269      ierr =        ice_alloc        ()      ! ice variables
270      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
271      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics
272      IF( ln_limdyn )   ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics
273      !
274      IF( lk_mpp    )   CALL mpp_sum( ierr )
275      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
276      !
277      CALL ice_rhg_init                ! set ice dynamics parameters
278      !
279      CALL ice_itd_init                ! ice thickness distribution initialization
280      !
281      CALL lim_thd_init                ! set ice thermodynics parameters
282      !
283      CALL lim_thd_sal_init            ! set ice salinity parameters
284       
285      ! MV MP 2016
286      CALL lim_mp_init                 ! set melt ponds parameters
287      ! END MV MP 2016
288
289      IF( ln_limdyn )   CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation
290      !                                ! Initial sea-ice state
291      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
292         CALL lim_istate
293      ELSE                                    ! start from a restart file
294         CALL lim_rst_read
295      ENDIF
296      CALL lim_var_agg(2)
297      CALL lim_var_glo2eqv
298      !
299      CALL lim_sbc_init                 ! ice surface boundary condition
300      !
301      IF( ln_limdiahsb) CALL lim_diahsb_init  ! initialization for diags
302      !
303      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction
304      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
305      !
306      DO jj = 1, jpj
307         DO ji = 1, jpi
308            IF( gphit(ji,jj) > 0._wp ) THEN  ;  rn_amax_2d(ji,jj) = rn_amax_n  ! NH
309            ELSE                             ;  rn_amax_2d(ji,jj) = rn_amax_s  ! SH
310            ENDIF
311         END DO
312      END DO
313      !
314   END SUBROUTINE ice_init
315
316
317   SUBROUTINE ice_run_init
318      !!-------------------------------------------------------------------
319      !!                  ***  ROUTINE ice_run_init ***
320      !!
321      !! ** Purpose :   Definition some run parameter for ice model
322      !!
323      !! ** Method  :   Read the namicerun namelist and check the parameter
324      !!              values called at the first timestep (nit000)
325      !!
326      !! ** input   :   Namelist namicerun
327      !!-------------------------------------------------------------------
328      INTEGER  ::   ios                 ! Local integer output status for namelist read
329      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   &
330         &                cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice 
331      NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt 
332      !!-------------------------------------------------------------------
333      !
334      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice
335      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
336901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
337
338      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice
339      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
340902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
341      IF(lwm) WRITE ( numoni, namicerun )
342      !
343      REWIND( numnam_ice_ref )              ! Namelist namicediag in reference namelist : Parameters for ice
344      READ  ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
345903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
346
347      REWIND( numnam_ice_cfg )              ! Namelist namicediag in configuration namelist : Parameters for ice
348      READ  ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
349904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
350      IF(lwm) WRITE ( numoni, namicediag )
351      !
352      IF(lwp) THEN                        ! control print
353         WRITE(numout,*)
354         WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
355         WRITE(numout,*) ' ~~~~~~'
356         WRITE(numout,*) '   number of ice  categories                              jpl    = ', jpl
357         WRITE(numout,*) '   number of ice  layers                                  nlay_i = ', nlay_i
358         WRITE(numout,*) '   number of snow layers                                  nlay_s = ', nlay_s
359         WRITE(numout,*) '   virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat
360         WRITE(numout,*) '   maximum ice concentration for NH                              = ', rn_amax_n 
361         WRITE(numout,*) '   maximum ice concentration for SH                              = ', rn_amax_s
362         WRITE(numout,*) '   Ice thermodynamics (T) or not (F)                  ln_limthd  = ', ln_limthd
363         WRITE(numout,*) '   Ice dynamics       (T) or not (F)                  ln_limdyn  = ', ln_limdyn
364         WRITE(numout,*) '     (ln_limdyn=T) Ice dynamics switch                nn_limdyn  = ', nn_limdyn
365         WRITE(numout,*) '       2: total'
366         WRITE(numout,*) '       1: advection only (no diffusion, no ridging/rafting)'
367         WRITE(numout,*) '       0: advection only (as 1 + prescribed velocity, bypass rheology)'
368         WRITE(numout,*) '     (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0)           = ', rn_uice
369         WRITE(numout,*) '     (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0)           = ', rn_vice
370         WRITE(numout,*)
371         WRITE(numout,*) '...and ice diagnostics'
372         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
373         WRITE(numout,*) '   Diagnose online heat/mass/salt budget     ln_limdiachk  = ', ln_limdiachk
374         WRITE(numout,*) '   Output          heat/mass/salt budget     ln_limdiahsb  = ', ln_limdiahsb
375         WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl
376         WRITE(numout,*) '   i-index for control prints (ln_limctl=true)             = ', iiceprt
377         WRITE(numout,*) '   j-index for control prints (ln_limctl=true)             = ', jiceprt
378      ENDIF
379      !
380      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
381         nn_monocat = 0
382         IF(lwp) WRITE(numout,*)
383         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
384      ENDIF
385      IF ( ( jpl == 1 ) .AND. ( nn_monocat == 0 ) ) THEN
386         CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
387      ENDIF
388      !
389      ! sea-ice timestep and inverse
390      rdt_ice   = REAL(nn_fsbc) * rdt 
391      r1_rdtice = 1._wp / rdt_ice 
392
393      ! inverse of nlay_i and nlay_s
394      r1_nlay_i = 1._wp / REAL( nlay_i, wp )
395      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
396      !
397      IF( lwp .AND. ln_bdy .AND. ln_limdiachk )  &
398         &   CALL ctl_warn('online conservation check activated but it does not work with BDY')
399      !
400      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice
401      !
402   END SUBROUTINE ice_run_init
403
404
405   SUBROUTINE ice_itd_init
406      !!------------------------------------------------------------------
407      !!                ***  ROUTINE ice_itd_init ***
408      !!
409      !! ** Purpose :   Initializes the ice thickness distribution
410      !! ** Method  :   ...
411      !! ** input   :   Namelist namiceitd
412      !!-------------------------------------------------------------------
413      INTEGER  ::   ios                 ! Local integer output status for namelist read
414      NAMELIST/namiceitd/ rn_himean
415      !
416      INTEGER  ::   jl                   ! dummy loop index
417      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
418      REAL(wp) ::   zhmax, znum, zden, zalpha !
419      !!------------------------------------------------------------------
420      !
421      REWIND( numnam_ice_ref )              ! Namelist namiceitd in reference namelist : Parameters for ice
422      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
423905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
424
425      REWIND( numnam_ice_cfg )              ! Namelist namiceitd in configuration namelist : Parameters for ice
426      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
427906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
428      IF(lwm) WRITE ( numoni, namiceitd )
429      !
430      IF(lwp) THEN                        ! control print
431         WRITE(numout,*)
432         WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution '
433         WRITE(numout,*) '~~~~~~~~~~~~'
434         WRITE(numout,*) '   mean ice thickness in the domain            rn_himean = ', rn_himean
435      ENDIF
436      !
437      !----------------------------------
438      !- Thickness categories boundaries
439      !----------------------------------
440      !
441      !==  h^(-alpha) function  ==!
442      zalpha = 0.05_wp
443      zhmax  = 3._wp * rn_himean
444      DO jl = 1, jpl
445         znum = jpl * ( zhmax+1 )**zalpha
446         zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
447         hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
448      END DO
449      !
450      DO jl = 1, jpl                ! mean thickness by category
451         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
452      END DO
453      !
454      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
455      !
456      IF(lwp) WRITE(numout,*)
457      IF(lwp) WRITE(numout,*) '      Thickness category boundaries '
458      IF(lwp) WRITE(numout,*) '         hi_max ', hi_max(0:jpl)
459      !
460   END SUBROUTINE ice_itd_init
461
462   SUBROUTINE ice_bef
463      !!----------------------------------------------------------------------
464      !!                  ***  ROUTINE ice_bef  ***
465      !!
466      !! ** purpose :  store ice variables at "before" time step
467      !!----------------------------------------------------------------------
468      INTEGER  ::   ji, jj, jl      ! dummy loop index
469      !!----------------------------------------------------------------------
470      !
471      DO jl = 1, jpl
472         DO jj = 2, jpjm1
473            DO ji = 2, jpim1
474               a_i_b  (ji,jj,jl)   = a_i  (ji,jj,jl)     ! ice area
475               v_i_b  (ji,jj,jl)   = v_i  (ji,jj,jl)     ! ice volume
476               v_s_b  (ji,jj,jl)   = v_s  (ji,jj,jl)     ! snow volume
477               smv_i_b(ji,jj,jl)   = smv_i(ji,jj,jl)     ! salt content
478               oa_i_b (ji,jj,jl)   = oa_i (ji,jj,jl)     ! areal age content
479               e_s_b  (ji,jj,:,jl) = e_s  (ji,jj,:,jl)   ! snow thermal energy
480               e_i_b  (ji,jj,:,jl) = e_i  (ji,jj,:,jl)   ! ice thermal energy
481               !                                         ! ice thickness
482               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes
483               ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
484               ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
485            END DO
486         END DO   
487      END DO
488      CALL lbc_lnk_multi( a_i_b, 'T', 1., v_i_b , 'T', 1., v_s_b , 'T', 1., smv_i_b, 'T', 1., &
489         &               oa_i_b, 'T', 1., ht_i_b, 'T', 1., ht_s_b, 'T', 1. )
490      CALL lbc_lnk( e_i_b, 'T', 1. )
491      CALL lbc_lnk( e_s_b, 'T', 1. )
492     
493      ! ice velocities & total concentration
494      DO jj = 2, jpjm1
495         DO ji = 2, jpim1
496            at_i_b(ji,jj)  = SUM( a_i_b(ji,jj,:) )
497            u_ice_b(ji,jj) = u_ice(ji,jj)
498            v_ice_b(ji,jj) = v_ice(ji,jj)
499         END DO
500      END DO
501      CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. )
502     
503   END SUBROUTINE ice_bef
504
505
506   SUBROUTINE ice_diag0
507      !!----------------------------------------------------------------------
508      !!                  ***  ROUTINE ice_diag0  ***
509      !!
510      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
511      !!               of the time step
512      !!----------------------------------------------------------------------
513      INTEGER  ::   ji, jj      ! dummy loop index
514      !!----------------------------------------------------------------------
515      DO jj = 1, jpj
516         DO ji = 1, jpi
517            sfx    (ji,jj) = 0._wp   ;
518            sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp
519            sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp
520            sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp
521            sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp
522            sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp
523            !
524            wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp
525            wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp
526            wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp
527            wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp
528            wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp
529            wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp 
530            wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
531            wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
532            wfx_snw_sni(ji,jj) = 0._wp 
533            ! MV MP 2016
534            wfx_pnd(ji,jj) = 0._wp
535            ! END MV MP 2016
536           
537            hfx_thd(ji,jj) = 0._wp   ;
538            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
539            hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp
540            hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp
541            hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp
542            hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp
543            hfx_err(ji,jj) = 0._wp   ;   hfx_err_rem(ji,jj) = 0._wp
544            hfx_err_dif(ji,jj) = 0._wp
545            wfx_err_sub(ji,jj) = 0._wp
546            !
547            afx_tot(ji,jj) = 0._wp   ;
548            afx_dyn(ji,jj) = 0._wp   ;   afx_thd(ji,jj) = 0._wp
549            !
550            diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp
551            diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp
552           
553            ! SIMIP diagnostics
554            diag_dms_dyn(ji,jj)  = 0._wp ; diag_dmi_dyn(ji,jj)  = 0._wp
555            diag_fc_bo(ji,jj)    = 0._wp ; diag_fc_su(ji,jj)    = 0._wp
556           
557            tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
558         END DO
559      END DO
560     
561   END SUBROUTINE ice_diag0
562
563#else
564   !!----------------------------------------------------------------------
565   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
566   !!----------------------------------------------------------------------
567CONTAINS
568   !
569   SUBROUTINE ice_stp ( kt, ksbc )     ! Dummy routine
570      INTEGER, INTENT(in) ::   kt, ksbc
571      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
572   END SUBROUTINE ice_stp
573   !
574   SUBROUTINE ice_init                 ! Dummy routine
575   END SUBROUTINE ice_init
576#endif
577
578   !!======================================================================
579END MODULE icestp
Note: See TracBrowser for help on using the repository browser.