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

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

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