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

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

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