source: NEMO/trunk/src/ICE/icestp.F90 @ 9604

Last change on this file since 9604 was 9604, checked in by clem, 2 years ago

change history of the ice routines

File size: 22.1 KB
Line 
1MODULE icestp
2   !!======================================================================
3   !!                       ***  MODULE  icestp  ***
4   !! sea ice : Master routine for all the sea ice model
5   !!=====================================================================
6   !!
7   !! The sea ice model SI3 (Sea Ice modelling Integrated Initiative),
8   !!                        aka Sea Ice cube for its nickname
9   !!
10   !!    is originally based on LIM3, developed in Louvain-la-Neuve by:
11   !!       * Martin Vancoppenolle (UCL-ASTR, Belgium)
12   !!       * Sylvain Bouillon (UCL-ASTR, Belgium)
13   !!       * Miguel Angel Morales Maqueda (NOC-L, UK)
14   !!      thanks to valuable earlier work by
15   !!       * Thierry Fichefet
16   !!       * Hugues Goosse
17   !!      thanks also to the following persons who contributed
18   !!       * Gurvan Madec, Claude Talandier, Christian Ethe (LOCEAN, France)
19   !!       * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany)
20   !!       * Bill Lipscomb (LANL), Cecilia Bitz (UWa) and Elisabeth Hunke (LANL), USA.
21   !!
22   !! SI3 has been made possible by a handful of persons who met as working group
23   !!     in France and UK:
24   !!    * ...
25   !!======================================================================
26   !! History :  4.0  !  2018     (C. Rousset)      Original code SI3
27   !!----------------------------------------------------------------------
28#if defined key_si3
29   !!----------------------------------------------------------------------
30   !!   'key_si3'                                       SI3 sea-ice model
31   !!----------------------------------------------------------------------
32   !!   ice_stp       : sea-ice model time-stepping and update ocean SBC over ice-covered area
33   !!   ice_init      : initialize sea-ice
34   !!----------------------------------------------------------------------
35   USE oce            ! ocean dynamics and tracers
36   USE dom_oce        ! ocean space and time domain
37   USE c1d            ! 1D vertical configuration
38   USE ice            ! sea-ice: variables
39   USE ice1D          ! sea-ice: thermodynamical 1D variables
40   !
41   USE phycst         ! Define parameters for the routines
42   USE eosbn2         ! equation of state
43   USE sbc_oce        ! Surface boundary condition: ocean fields
44   USE sbc_ice        ! Surface boundary condition: ice   fields
45   !
46   USE iceforcing     ! sea-ice: Surface boundary condition       !!gm why not icesbc module name
47   USE icedyn         ! sea-ice: dynamics
48   USE icethd         ! sea-ice: thermodynamics
49   USE icecor         ! sea-ice: corrections
50   USE iceupdate      ! sea-ice: sea surface boundary condition update
51   USE icedia         ! sea-ice: budget diagnostics
52   USE icewri         ! sea-ice: outputs
53   USE icerst         ! sea-ice: restarts
54   USE icevar         ! sea-ice: operations
55   USE icectl         ! sea-ice: control
56   USE iceistate      ! sea-ice: initial state
57   USE iceitd         ! sea-ice: remapping thickness distribution
58   USE icealb         ! sea-ice: albedo
59   !
60   USE bdy_oce , ONLY : ln_bdy   ! flag for bdy
61   USE bdyice         ! unstructured open boundary data for sea-ice
62# if defined key_agrif
63   USE agrif_ice
64   USE agrif_ice_interp
65# endif
66   !
67   USE in_out_manager ! I/O manager
68   USE iom            ! I/O manager library
69   USE lib_mpp        ! MPP library
70   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
71   USE timing         ! Timing
72   USE prtctl         ! Print control
73
74   IMPLICIT NONE
75   PRIVATE
76
77   PUBLIC   ice_stp    ! called by sbcmod.F90
78   PUBLIC   ice_init   ! called by sbcmod.F90
79
80   !! * Substitutions
81#  include "vectopt_loop_substitute.h90"
82   !!----------------------------------------------------------------------
83   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
84   !! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $
85   !! Software governed by the CeCILL licence     (./LICENSE)
86   !!----------------------------------------------------------------------
87CONTAINS
88
89   SUBROUTINE ice_stp( kt, ksbc )
90      !!---------------------------------------------------------------------
91      !!                  ***  ROUTINE ice_stp  ***
92      !!
93      !! ** Purpose :   sea-ice model time-stepping and update ocean surface
94      !!              boundary condition over ice-covered area
95      !!
96      !! ** Method  :   ice model time stepping
97      !!              - call the ice dynamics routine
98      !!              - call the ice advection/diffusion routine
99      !!              - call the ice thermodynamics routine
100      !!              - call the routine that computes mass and
101      !!                heat fluxes at the ice/ocean interface
102      !!              - save the outputs
103      !!              - save the outputs for restart when necessary
104      !!
105      !! ** Action  : - time evolution of the LIM sea-ice model
106      !!              - update all sbc variables below sea-ice:
107      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx
108      !!---------------------------------------------------------------------
109      INTEGER, INTENT(in) ::   kt      ! ocean time step
110      INTEGER, INTENT(in) ::   ksbc    ! flux formulation (user defined, bulk, or Pure Coupled)
111      !
112      INTEGER ::   jl   ! dummy loop index
113      !!----------------------------------------------------------------------
114      !
115      IF( ln_timing )   CALL timing_start('ice_stp')
116      !
117      !                                      !-----------------------!
118      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN   ! --- Ice time step --- !
119         !                                   !-----------------------!
120         !
121         kt_ice = kt                              ! -- Ice model time step
122         !
123         u_oce(:,:) = ssu_m(:,:)                  ! -- mean surface ocean current
124         v_oce(:,:) = ssv_m(:,:)
125         !
126         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )   ! -- freezing temperature [Kelvin] (set to rt0 over land)
127         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
128         !
129         !                          !==  AGRIF Parent to Child  ==!
130#if defined key_agrif
131         !                              ! lim_nbstep ranges from 1 to the nb of child ocean steps inside one parent ice step
132         IF( .NOT. Agrif_Root() )       lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
133         !                              ! these calls must remain here for restartability purposes
134                                        CALL agrif_interp_si3( 'T' ) 
135                                        CALL agrif_interp_si3( 'U' )
136                                        CALL agrif_interp_si3( 'V' )
137#endif
138                                        CALL store_fields             ! Store now ice values
139         !
140         !------------------------------------------------!
141         ! --- Dynamical coupling with the atmosphere --- !
142         !------------------------------------------------!
143         ! It provides the following fields used in sea ice model:
144         !    utau_ice, vtau_ice = surface ice stress [N/m2]
145         !------------------------------------------------!
146                                        CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )         
147         !-------------------------------------!
148         ! --- ice dynamics and advection  --- !
149         !-------------------------------------!
150                                        CALL diag_set0                ! set diag of mass, heat and salt fluxes to 0
151                                        CALL ice_rst_opn( kt )        ! Open Ice restart file (if necessary)
152         !
153         IF( ln_icedyn .AND. .NOT.lk_c1d )   &
154            &                           CALL ice_dyn( kt )            ! -- Ice dynamics
155         !
156         !                          !==  lateral boundary conditions  ==!
157         IF( ln_icethd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
158         !
159         !                          !==  previous lead fraction and ice volume for flux calculations
160                                        CALL ice_var_glo2eqv          ! h_i and h_s for ice albedo calculation
161                                        CALL ice_var_agg(1)           ! at_i for coupling
162                                        CALL store_fields             ! Store now ice values
163         !
164         !------------------------------------------------------!
165         ! --- Thermodynamical coupling with the atmosphere --- !
166         !------------------------------------------------------!
167         ! It provides the following fields used in sea ice model:
168         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s]
169         !    sprecip              = solid precipitation                           [Kg/m2/s]
170         !    evap_ice             = sublimation                                   [Kg/m2/s]
171         !    qsr_tot , qns_tot    = solar & non solar heat flux (total)           [W/m2]
172         !    qsr_ice , qns_ice    = solar & non solar heat flux over ice          [W/m2]
173         !    dqns_ice             = non solar  heat sensistivity                  [W/m2]
174         !    qemp_oce, qemp_ice,  = sensible heat (associated with evap & precip) [W/m2]
175         !    qprec_ice, qevap_ice
176         !------------------------------------------------------!
177                                        CALL ice_forcing_flx( kt, ksbc )
178         !----------------------------!
179         ! --- ice thermodynamics --- !
180         !----------------------------!
181         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics     
182         !
183         !
184         IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections
185         !
186                                        CALL ice_var_glo2eqv          ! necessary calls (at least for coupling)
187                                        CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling)
188         !
189!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
190!!       moreover it should only be called at the update frequency (as in agrif_ice_update.F90)
191!# if defined key_agrif
192!         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid()
193!# endif
194                                        CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
195!# if defined key_agrif
196!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid()
197!# endif
198         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs
199         !
200                                        CALL ice_wri( kt )            ! -- Ice outputs
201         !
202         IF( lrst_ice )                 CALL ice_rst_write( kt )      ! -- Ice restart file
203         !
204         IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash
205         !
206      ENDIF   ! End sea-ice time step only
207
208      !-------------------------!
209      ! --- Ocean time step --- !
210      !-------------------------!
211      IF( ln_icedyn )                   CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses
212!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
213      !
214      IF( ln_timing )   CALL timing_stop('ice_stp')
215      !
216   END SUBROUTINE ice_stp
217
218
219   SUBROUTINE ice_init
220      !!----------------------------------------------------------------------
221      !!                  ***  ROUTINE ice_init  ***
222      !!
223      !! ** purpose :   Initialize sea-ice parameters
224      !!----------------------------------------------------------------------
225      INTEGER :: ji, jj, ierr
226      !!----------------------------------------------------------------------
227      IF(lwp) WRITE(numout,*)
228      IF(lwp) WRITE(numout,*) 'ice_init: Arrays allocation & Initialization off all routines & init state' 
229      IF(lwp) WRITE(numout,*) '~~~~~~~~'
230      !
231      !                                ! Open the reference and configuration namelist files and namelist output file
232      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
233      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
234      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
235      !
236      CALL par_init                ! set some ice run parameters
237      !
238      !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
239      ierr =        ice_alloc        ()      ! ice variables
240      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
241      ierr = ierr + ice1D_alloc      ()      ! thermodynamics
242      !
243      IF( lk_mpp    )   CALL mpp_sum( ierr )
244      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
245      !
246      CALL ice_itd_init                ! ice thickness distribution initialization
247      !
248      CALL ice_thd_init                ! set ice thermodynics parameters (clem: important to call it first for melt ponds)
249      !
250      !                                ! Initial sea-ice state
251      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
252         CALL ice_istate_init
253         CALL ice_istate
254      ELSE                                    ! start from a restart file
255         CALL ice_rst_read
256      ENDIF
257      CALL ice_var_glo2eqv
258      CALL ice_var_agg(2)
259      !
260      CALL ice_forcing_init            ! set ice-ocean and ice-atm. coupling parameters
261      !
262      CALL ice_dyn_init                ! set ice dynamics parameters
263      !
264      CALL ice_update_init             ! ice surface boundary condition
265      !
266      CALL ice_alb_init                ! ice surface albedo
267      !
268      CALL ice_dia_init                ! initialization for diags
269      !
270      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction
271      tn_ice(:,:,:) = t_su(:,:,:)      ! initialisation of surface temp for coupled simu
272      !
273      !                                ! set max concentration in both hemispheres
274      WHERE( gphit(:,:) > 0._wp )   ;   rn_amax_2d(:,:) = rn_amax_n  ! NH
275      ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH
276      END WHERE
277
278      IF( ln_rstart )   CALL iom_close( numrir )  ! close input ice restart file
279      !
280   END SUBROUTINE ice_init
281
282
283   SUBROUTINE par_init
284      !!-------------------------------------------------------------------
285      !!                  ***  ROUTINE par_init ***
286      !!
287      !! ** Purpose :   Definition generic parameters for ice model
288      !!
289      !! ** Method  :   Read namelist and check the parameter
290      !!                values called at the first timestep (nit000)
291      !!
292      !! ** input   :   Namelist nampar
293      !!-------------------------------------------------------------------
294      INTEGER  ::   ios                 ! Local integer
295      !!
296      NAMELIST/nampar/ jpl, nlay_i, nlay_s, nn_virtual_itd, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s,  &
297         &             cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir
298      !!-------------------------------------------------------------------
299      !
300      REWIND( numnam_ice_ref )      ! Namelist nampar in reference namelist : Parameters for ice
301      READ  ( numnam_ice_ref, nampar, IOSTAT = ios, ERR = 901)
302901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampar in reference namelist', lwp )
303      REWIND( numnam_ice_cfg )      ! Namelist nampar in configuration namelist : Parameters for ice
304      READ  ( numnam_ice_cfg, nampar, IOSTAT = ios, ERR = 902 )
305902   IF( ios > 0 )   CALL ctl_nam ( ios , 'nampar in configuration namelist', lwp )
306      IF(lwm) WRITE( numoni, nampar )
307      !
308      IF(lwp) THEN                  ! control print
309         WRITE(numout,*)
310         WRITE(numout,*) '   par_init: ice parameters shared among all the routines'
311         WRITE(numout,*) '   ~~~~~~~~'
312         WRITE(numout,*) '      Namelist nampar: '
313         WRITE(numout,*) '         number of ice  categories                           jpl       = ', jpl
314         WRITE(numout,*) '         number of ice  layers                               nlay_i    = ', nlay_i
315         WRITE(numout,*) '         number of snow layers                               nlay_s    = ', nlay_s
316         WRITE(numout,*) '         virtual ITD param for jpl=1 (1-3) or not (0)   nn_virtual_itd = ', nn_virtual_itd
317         WRITE(numout,*) '         Ice dynamics       (T) or not (F)                   ln_icedyn = ', ln_icedyn
318         WRITE(numout,*) '         Ice thermodynamics (T) or not (F)                   ln_icethd = ', ln_icethd
319         WRITE(numout,*) '         maximum ice concentration for NH                              = ', rn_amax_n 
320         WRITE(numout,*) '         maximum ice concentration for SH                              = ', rn_amax_s
321      ENDIF
322      !                                        !--- check consistency
323      IF ( jpl > 1 .AND. nn_virtual_itd == 1 ) THEN
324         nn_virtual_itd = 0
325         IF(lwp) WRITE(numout,*)
326         IF(lwp) WRITE(numout,*) '   nn_virtual_itd forced to 0 as jpl>1, no need with multiple categories to emulate them'
327      ENDIF
328      !
329      IF( ln_cpl .AND. nn_cats_cpl /= 1 .AND. nn_cats_cpl /= jpl ) THEN
330         CALL ctl_stop( 'STOP', 'par_init: in coupled mode, nn_cats_cpl should be either 1 or jpl' )
331      ENDIF
332      !
333      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY')
334      !
335      rdt_ice   = REAL(nn_fsbc) * rdt          !--- sea-ice timestep and its inverse
336      r1_rdtice = 1._wp / rdt_ice
337      IF(lwp) WRITE(numout,*)
338      IF(lwp) WRITE(numout,*) '      ice timestep rdt_ice = nn_fsbc*rdt = ', rdt_ice
339      !
340      r1_nlay_i = 1._wp / REAL( nlay_i, wp )   !--- inverse of nlay_i and nlay_s
341      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
342      !
343   END SUBROUTINE par_init
344
345
346   SUBROUTINE store_fields
347      !!----------------------------------------------------------------------
348      !!                  ***  ROUTINE store_fields  ***
349      !!
350      !! ** purpose :  store ice variables at "before" time step
351      !!----------------------------------------------------------------------
352      INTEGER  ::   ji, jj, jl      ! dummy loop index
353      !!----------------------------------------------------------------------
354      !
355      a_i_b (:,:,:)   = a_i (:,:,:)     ! ice area
356      v_i_b (:,:,:)   = v_i (:,:,:)     ! ice volume
357      v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume
358      sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content
359      oa_i_b(:,:,:)   = oa_i(:,:,:)     ! areal age content
360      e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy
361      e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy
362      WHERE( a_i_b(:,:,:) >= epsi20 )
363         h_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:)   ! ice thickness
364         h_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:)   ! snw thickness
365      ELSEWHERE
366         h_i_b(:,:,:) = 0._wp
367         h_s_b(:,:,:) = 0._wp
368      END WHERE
369      !
370      ! ice velocities & total concentration
371      at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 )
372      u_ice_b(:,:) = u_ice(:,:)
373      v_ice_b(:,:) = v_ice(:,:)
374      !
375   END SUBROUTINE store_fields
376
377
378   SUBROUTINE diag_set0
379      !!----------------------------------------------------------------------
380      !!                  ***  ROUTINE diag_set0  ***
381      !!
382      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
383      !!               of the time step
384      !!----------------------------------------------------------------------
385      INTEGER  ::   ji, jj      ! dummy loop index
386      !!----------------------------------------------------------------------
387      sfx    (:,:) = 0._wp   ;
388      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp
389      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
390      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
391      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
392      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
393      !
394      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
395      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
396      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
397      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
398      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
399      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp 
400      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
401      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
402      wfx_snw_sni(:,:) = 0._wp 
403      wfx_pnd(:,:) = 0._wp
404
405      hfx_thd(:,:) = 0._wp   ;
406      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
407      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
408      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
409      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
410      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
411      hfx_err_rem(:,:) = 0._wp
412      hfx_err_dif(:,:) = 0._wp
413      wfx_err_sub(:,:) = 0._wp
414      !
415      afx_tot(:,:) = 0._wp   ;
416      !
417      diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp
418      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
419
420      ! SIMIP diagnostics
421      diag_fc_bo(:,:) = 0._wp ; diag_fc_su(:,:) = 0._wp
422      t_si(:,:,:) = rt0       ! temp at the ice-snow interface
423
424      tau_icebfr(:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here)
425      cnd_ice   (:,:,:) = 0._wp   ! initialisation of the effective conductivity at the top of ice/snow (Jules coupling)
426      !
427      ! for control checks (ln_icediachk)
428      diag_trp_vi(:,:) = 0._wp   ;   diag_trp_vs(:,:) = 0._wp
429      diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp
430      diag_trp_sv(:,:) = 0._wp
431
432   END SUBROUTINE diag_set0
433
434#else
435   !!----------------------------------------------------------------------
436   !!   Default option           Dummy module         NO SI3 sea-ice model
437   !!----------------------------------------------------------------------
438CONTAINS
439   SUBROUTINE ice_stp ( kt, ksbc )     ! Dummy routine
440      INTEGER, INTENT(in) ::   kt, ksbc
441      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
442   END SUBROUTINE ice_stp
443   SUBROUTINE ice_init                 ! Dummy routine
444   END SUBROUTINE ice_init
445#endif
446
447   !!======================================================================
448END MODULE icestp
Note: See TracBrowser for help on using the repository browser.