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

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

changes in style - part1 - (now the code looks better txs to Gurvan's comments)
File size: 30.1 KB
Line 
1MODULE icestp
2   !!======================================================================
3   !!                       ***  MODULE  icestp  ***
4   !! Surface module :  update the ocean surface boundary condition over ice
5   !!       &           covered area using LIM sea-ice model
6   !! Sea-Ice model  :  LIM-3 Sea ice model time-stepping
7   !!=====================================================================
8   !! History :  2.0  ! 2006-12  (M. Vancoppenolle) Original code
9   !!            3.0  ! 2008-02  (C. Talandier)  Surface module from icestp.F90
10   !!             -   ! 2008-04  (G. Madec)  sltyle and 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
130         !                 ! mean surface ocean current at ice velocity point
131         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
132         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
133!!gm the umask, vmask above is useless as ssu_m, ssv_m are build from masked un,vn (used t be here for B-grid sea-ice)
134
135         !                 ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
136         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
137         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
138         !
139         CALL ice_bef      ! Store previous ice values
140         !
141         !------------------------------------------------!
142         ! --- Dynamical coupling with the atmosphere --- !
143         !------------------------------------------------!
144         ! It provides the following fields used in sea ice model:
145         !    utau_ice, vtau_ice = surface ice stress [N/m2]
146         !------------------------------------------------!
147         CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )
148                                     
149         !-------------------------------------!
150         ! --- ice dynamics and advection  --- !
151         !-------------------------------------!
152         CALL ice_diag0             ! set diag of mass, heat and salt fluxes to 0
153         CALL ice_rst_opn( kt )     ! Open Ice restart file (if necessary)
154         !
155         SELECT CASE( nice_dyn )   
156         CASE ( np_dynFULL )          !==  all dynamical processes  ==!
157            CALL ice_rhg   ( kt )          ! -- rheology 
158            CALL ice_adv   ( kt )          ! -- advection of ice
159            CALL ice_rdgrft( kt )          ! -- ridging/rafting
160            CALL ice_cor   ( kt , 1 )      ! -- Corrections
161         CASE ( np_dyn )              !==  pure dynamics only ==!   (no ridging/rafting)   (nono cat. case 2)
162            CALL ice_rhg   ( kt )          ! -- rheology 
163            CALL ice_adv   ( kt )          ! -- advection of ice
164            CALL ice_cor   ( kt , 1 )      ! -- Corrections
165         CASE ( np_dynPURE )          !==  pure dynamics only ==!   (nn_limdyn= 0 or 1 )
166            CALL ice_rhg   ( kt )          ! -- rheology 
167            CALL ice_adv   ( kt )          ! -- advection of ice
168         END SELECT
169         !
170         !                          !==  lateral boundary conditions  ==!
171#if defined key_agrif
172         IF( .NOT. Agrif_Root()     )   CALL agrif_interp_lim3('T')   ! -- AGRIF
173#endif
174         IF( ln_limthd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
175         !
176         !
177         !                          !==  previous lead fraction and ice volume for flux calculations
178         !
179         CALL ice_var_glo2eqv            ! ht_i and ht_s for ice albedo calculation
180         CALL ice_var_agg(1)             ! at_i for coupling
181         CALL ice_bef                    ! Store previous ice values
182
183         !------------------------------------------------------!
184         ! --- Thermodynamical coupling with the atmosphere --- !
185         !------------------------------------------------------!
186         ! It provides the following fields used in sea ice model:
187         !    fr1_i0  , fr2_i0     = 1sr & 2nd fraction of qsr penetration in ice  [%]
188         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s]
189         !    sprecip              = solid precipitation                           [Kg/m2/s]
190         !    evap_ice             = sublimation                                   [Kg/m2/s]
191         !    qsr_tot , qns_tot    = solar & non solar heat flux (total)           [W/m2]
192         !    qsr_ice , qns_ice    = solar & non solar heat flux over ice          [W/m2]
193         !    dqns_ice             = non solar  heat sensistivity                  [W/m2]
194         !    qemp_oce, qemp_ice,  = sensible heat (associated with evap & precip) [W/m2]
195         !    qprec_ice, qevap_ice
196         !------------------------------------------------------!
197         CALL ice_forcing_flx( kt, ksbc )
198
199         !----------------------------!
200         ! --- ice thermodynamics --- !
201         !----------------------------!
202         IF( ln_limthd )            CALL ice_thd( kt )          ! -- Ice thermodynamics     
203
204         ! MV MP 2016
205         IF ( ln_pnd )              CALL lim_mp( kt )           ! -- Melt ponds
206         ! END MV MP 2016
207
208         IF( ln_limthd )            CALL ice_cor( kt , 2 )      ! -- Corrections
209         ! ---
210# if defined key_agrif
211         IF( .NOT. Agrif_Root() )   CALL agrif_update_lim3( kt )
212# endif
213                                    CALL ice_var_glo2eqv        ! necessary calls (at least for coupling)
214                                    CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling)
215                                    !
216!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
217!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
218!!# if defined key_agrif
219!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid()
220!!# endif
221                                    CALL ice_update_flx( kt )   ! -- Update ocean surface mass, heat and salt fluxes
222!!# if defined key_agrif
223!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
224!!# endif
225         IF( ln_limdiahsb )           CALL ice_dia( kt )        ! -- Diagnostics and outputs
226         !
227                                      CALL ice_wri( 1 )         ! -- Ice outputs
228         !
229         IF( kt == nit000 .AND. ln_rstart )   &                !!gm  I don't understand the ln_rstart, if needed, add a comment, please
230            &                         CALL iom_close( numrir )  ! close input ice restart file
231         !
232         IF( lrst_ice )               CALL ice_rst_write( kt )  ! -- Ice restart file
233         !
234         IF( ln_limctl )              CALL ice_ctl( kt )        ! alerts in case of model crash
235         !
236      ENDIF   ! End sea-ice time step only
237
238      !-------------------------!
239      ! --- Ocean time step --- !
240      !-------------------------!
241      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
242      !    using before instantaneous surf. currents
243      IF( ln_limdyn )                 CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )
244!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
245      !
246      IF( nn_timing == 1 ) CALL timing_stop('ice_stp')
247      !
248   END SUBROUTINE ice_stp
249
250
251   SUBROUTINE ice_init
252      !!----------------------------------------------------------------------
253      !!                  ***  ROUTINE ice_init  ***
254      !!
255      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules
256      !!----------------------------------------------------------------------
257      INTEGER :: ji, jj, ierr
258      !!----------------------------------------------------------------------
259      IF(lwp) WRITE(numout,*)
260      IF(lwp) WRITE(numout,*) 'ice_init : update ocean surface boudary condition' 
261      IF(lwp) WRITE(numout,*) '~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
262      !
263      !                                ! Open the reference and configuration namelist files and namelist output file
264      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
265      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
266      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
267      !
268      CALL ice_run_init                ! set some ice run parameters
269      !
270      !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
271      ierr =        ice_alloc        ()      ! ice variables
272      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
273      ierr = ierr + ice1D_alloc      ()      ! thermodynamics
274      IF( ln_limdyn )   ierr = ierr + ice_rdgrft_alloc ()      ! ridging/rafting
275      !
276      IF( lk_mpp    )   CALL mpp_sum( ierr )
277      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
278      !
279      CALL ice_rhg_init                ! set ice dynamics parameters
280      !
281      CALL ice_itd_init                ! ice thickness distribution initialization
282      !
283      CALL ice_thd_init                ! set ice thermodynics parameters
284      !
285      CALL ice_thd_sal_init            ! set ice salinity parameters
286       
287      ! MV MP 2016
288      CALL lim_mp_init                 ! set melt ponds parameters
289      ! END MV MP 2016
290
291      IF( ln_limdyn )   CALL ice_rdgrft_init             ! ice thickness distribution initialization for ridging/rafting
292      !                                ! Initial sea-ice state
293      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
294         CALL ice_ist
295      ELSE                                    ! start from a restart file
296         CALL ice_rst_read
297      ENDIF
298      CALL ice_var_agg(2)
299      CALL ice_var_glo2eqv
300      !
301      CALL ice_update_init                 ! ice surface boundary condition
302      !
303      IF( ln_limdiahsb) CALL ice_dia_init  ! initialization for diags
304      !
305      fr_i  (:,:)   = at_i(:,:)         ! initialisation of sea-ice fraction
306      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
307      !
308      DO jj = 1, jpj
309         DO ji = 1, jpi
310            IF( gphit(ji,jj) > 0._wp ) THEN   ;   rn_amax_2d(ji,jj) = rn_amax_n  ! NH
311            ELSE                              ;   rn_amax_2d(ji,jj) = rn_amax_s  ! SH
312            ENDIF
313         END DO
314      END DO
315      !
316   END SUBROUTINE ice_init
317
318
319   SUBROUTINE ice_run_init
320      !!-------------------------------------------------------------------
321      !!                  ***  ROUTINE ice_run_init ***
322      !!
323      !! ** Purpose :   Definition some run parameter for ice model
324      !!
325      !! ** Method  :   Read the namicerun namelist and check the parameter
326      !!              values called at the first timestep (nit000)
327      !!
328      !! ** input   :   Namelist namicerun
329      !!-------------------------------------------------------------------
330      INTEGER  ::   ios                 ! Local integer output status for namelist read
331      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   &
332         &                cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice 
333      NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt 
334      !!-------------------------------------------------------------------
335      !
336      REWIND( numnam_ice_ref )      ! Namelist namicerun in reference namelist : Parameters for ice
337      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
338901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
339
340      REWIND( numnam_ice_cfg )      ! Namelist namicerun in configuration namelist : Parameters for ice
341      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
342902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
343      IF(lwm) WRITE ( numoni, namicerun )
344      !
345      REWIND( numnam_ice_ref )      ! Namelist namicediag in reference namelist : Parameters for ice
346      READ  ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
347903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
348
349      REWIND( numnam_ice_cfg )      ! Namelist namicediag in configuration namelist : Parameters for ice
350      READ  ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
351904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
352      IF(lwm) WRITE ( numoni, namicediag )
353      !
354      IF(lwp) THEN                  ! control print
355         WRITE(numout,*)
356         WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
357         WRITE(numout,*) ' ~~~~~~'
358         WRITE(numout,*) '   Namelist namicerun : '
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,*) '         associated switch                               nn_limdyn  = ', nn_limdyn
368         WRITE(numout,*) '            =2 all processes (default option)'
369         WRITE(numout,*) '            =1 advection only (no ridging/rafting)'
370         WRITE(numout,*) '            =0 advection only with prescribed velocity given by '
371         WRITE(numout,*) '               a uniform field        (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')'
372         WRITE(numout,*)
373         WRITE(numout,*) '   Namelist namicediag : '
374         WRITE(numout,*) '      Diagnose online heat/mass/salt budget      ln_limdiachk = ', ln_limdiachk
375         WRITE(numout,*) '      Output          heat/mass/salt budget      ln_limdiahsb = ', ln_limdiahsb
376         WRITE(numout,*) '      control prints for a given grid point         ln_limctl = ', ln_limctl
377         WRITE(numout,*) '         chosen grid point position         (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')'
378      ENDIF
379      !
380      !                                         !--- check consistency
381      IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN
382         nn_monocat = 0
383         IF(lwp) WRITE(numout,*)
384         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
385      ENDIF
386      IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
387         CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
388      ENDIF
389      !
390      IF( ln_bdy .AND. ln_limdiachk )   CALL ctl_warn('ice_run_init: online conservation check does not work with BDY')
391      !
392      !                             ! set the choice of ice dynamics
393      IF( lk_c1d .OR. .NOT.ln_limdyn ) THEN
394         nice_dyn = np_dynNO                    !--- no dynamics
395      ELSE
396         SELECT CASE( nn_limdyn )
397         CASE( 2 )                   
398            IF( nn_monocat /= 2 ) THEN          !--- full dynamics (rheology + advection + ridging/rafting + correction)
399               nice_dyn = np_dynFULL
400            ELSE
401               nice_dyn = np_dyn                !--- dynamics without ridging/rafting
402            ENDIF
403         CASE( 0 , 1 )                          !--- dynamics without ridging/rafting and correction
404            nice_dyn = np_dynPURE
405         END SELECT
406      ENDIF
407      !                                         !--- simple conservative piling, comparable with LIM2
408      l_piling = nn_limdyn == 1 .OR. ( nn_monocat == 2  .AND.  jpl == 1 )
409      !
410      rdt_ice   = REAL(nn_fsbc) * rdt           !--- sea-ice timestep and inverse
411      r1_rdtice = 1._wp / rdt_ice
412      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice
413      !
414      r1_nlay_i = 1._wp / REAL( nlay_i, wp )    !--- inverse of nlay_i and nlay_s
415      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
416      !
417   END SUBROUTINE ice_run_init
418
419
420   SUBROUTINE ice_itd_init
421      !!------------------------------------------------------------------
422      !!                ***  ROUTINE ice_itd_init ***
423      !!
424      !! ** Purpose :   Initializes the ice thickness distribution
425      !! ** Method  :   ...
426      !! ** input   :   Namelist namiceitd
427      !!-------------------------------------------------------------------
428      INTEGER  ::   jl    ! dummy loop index
429      INTEGER  ::   ios   ! Local integer output status for namelist read
430      REAL(wp) ::   zc1, zc2, zc3, zx1          ! local scalars
431      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
432      !!
433      NAMELIST/namiceitd/ rn_himean
434      !!------------------------------------------------------------------
435      !
436      REWIND( numnam_ice_ref )      ! Namelist namiceitd in reference namelist : Parameters for ice
437      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
438905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
439
440      REWIND( numnam_ice_cfg )      ! Namelist namiceitd in configuration namelist : Parameters for ice
441      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
442906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
443      IF(lwm) WRITE ( numoni, namiceitd )
444      !
445      IF(lwp) THEN                  ! control print
446         WRITE(numout,*)
447         WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution '
448         WRITE(numout,*) '~~~~~~~~~~~~'
449         WRITE(numout,*) '   Namelist namicerun : '
450         WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean
451      ENDIF
452      !
453      !-----------------------------------!
454      !  Thickness categories boundaries  !
455      !-----------------------------------!
456      !
457      zalpha = 0.05_wp              ! max of each category (from h^(-alpha) function)
458      zhmax  = 3._wp * rn_himean
459      DO jl = 1, jpl
460         znum = jpl * ( zhmax+1 )**zalpha
461         zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
462         hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
463      END DO
464      !
465      DO jl = 1, jpl                ! mean thickness by category
466         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
467      END DO
468      !
469      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
470      !
471      IF(lwp) WRITE(numout,*)
472      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :'
473      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl)
474      !
475   END SUBROUTINE ice_itd_init
476
477
478   SUBROUTINE ice_bef
479      !!----------------------------------------------------------------------
480      !!                  ***  ROUTINE ice_bef  ***
481      !!
482      !! ** purpose :  store ice variables at "before" time step
483      !!----------------------------------------------------------------------
484      INTEGER  ::   ji, jj, jl      ! dummy loop index
485      !!----------------------------------------------------------------------
486      !
487      DO jl = 1, jpl
488         DO jj = 2, jpjm1
489            DO ji = 2, jpim1
490               a_i_b  (ji,jj,jl)   = a_i  (ji,jj,jl)     ! ice area
491               v_i_b  (ji,jj,jl)   = v_i  (ji,jj,jl)     ! ice volume
492               v_s_b  (ji,jj,jl)   = v_s  (ji,jj,jl)     ! snow volume
493               smv_i_b(ji,jj,jl)   = smv_i(ji,jj,jl)     ! salt content
494               oa_i_b (ji,jj,jl)   = oa_i (ji,jj,jl)     ! areal age content
495               e_s_b  (ji,jj,:,jl) = e_s  (ji,jj,:,jl)   ! snow thermal energy
496               e_i_b  (ji,jj,:,jl) = e_i  (ji,jj,:,jl)   ! ice thermal energy
497               !                                         ! ice thickness
498               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes
499               ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
500               ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
501            END DO
502         END DO   
503      END DO
504      CALL lbc_lnk_multi(  a_i_b, 'T', 1., v_i_b , 'T', 1., ht_i_b, 'T', 1., smv_i_b, 'T', 1.,   &
505         &                oa_i_b, 'T', 1., v_s_b , 'T', 1., ht_s_b, 'T', 1. )
506      CALL lbc_lnk( e_i_b, 'T', 1. )
507      CALL lbc_lnk( e_s_b, 'T', 1. )
508     
509!!gm Question:  here , a_i_b, u_ice and v_ice  are defined over the whole domain,
510!!gm            so why not just a copy over the whole domain and no lbc_lnk ?
511!!gm            that is some thing like:
512!            at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 )
513!            u_ice_b(:,:) = u_ice(:,:)
514!            v_ice_b(:,:) = v_ice(:,:)
515!    idem for the loop above
516!!gm
517      ! ice velocities & total concentration
518      DO jj = 2, jpjm1
519         DO ji = 2, jpim1
520            at_i_b(ji,jj)  = SUM( a_i_b(ji,jj,:) )
521            u_ice_b(ji,jj) = u_ice(ji,jj)
522            v_ice_b(ji,jj) = v_ice(ji,jj)
523         END DO
524      END DO
525      CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. )
526      !
527   END SUBROUTINE ice_bef
528
529
530   SUBROUTINE ice_diag0
531      !!----------------------------------------------------------------------
532      !!                  ***  ROUTINE ice_diag0  ***
533      !!
534      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
535      !!               of the time step
536      !!----------------------------------------------------------------------
537      INTEGER  ::   ji, jj      ! dummy loop index
538      !!----------------------------------------------------------------------
539      DO jj = 1, jpj
540         DO ji = 1, jpi
541            sfx    (ji,jj) = 0._wp   ;
542            sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp
543            sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp
544            sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp
545            sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp
546            sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp
547            !
548            wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp
549            wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp
550            wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp
551            wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp
552            wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp
553            wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp 
554            wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
555            wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
556            wfx_snw_sni(ji,jj) = 0._wp 
557            ! MV MP 2016
558            wfx_pnd(ji,jj) = 0._wp
559            ! END MV MP 2016
560           
561            hfx_thd(ji,jj) = 0._wp   ;
562            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
563            hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp
564            hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp
565            hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp
566            hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp
567            hfx_err(ji,jj) = 0._wp   ;   hfx_err_rem(ji,jj) = 0._wp
568            hfx_err_dif(ji,jj) = 0._wp
569            wfx_err_sub(ji,jj) = 0._wp
570            !
571            afx_tot(ji,jj) = 0._wp   ;
572            afx_dyn(ji,jj) = 0._wp   ;   afx_thd(ji,jj) = 0._wp
573            !
574            diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp
575            diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp
576           
577            ! SIMIP diagnostics
578            diag_dms_dyn(ji,jj)  = 0._wp ; diag_dmi_dyn(ji,jj)  = 0._wp
579            diag_fc_bo(ji,jj)    = 0._wp ; diag_fc_su(ji,jj)    = 0._wp
580           
581            tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
582         END DO
583      END DO
584     
585   END SUBROUTINE ice_diag0
586
587#else
588   !!----------------------------------------------------------------------
589   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
590   !!----------------------------------------------------------------------
591CONTAINS
592   SUBROUTINE ice_stp ( kt, ksbc )     ! Dummy routine
593      INTEGER, INTENT(in) ::   kt, ksbc
594      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
595   END SUBROUTINE ice_stp
596   SUBROUTINE ice_init                 ! Dummy routine
597   END SUBROUTINE ice_init
598#endif
599
600   !!======================================================================
601END MODULE icestp
Note: See TracBrowser for help on using the repository browser.