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_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90 @ 9012

Last change on this file since 9012 was 8933, checked in by clem, 6 years ago

dev_CNRS_2017: raw commit of the Met-Office specifics for coupling with their Earth system model. It is not yet tested and debuging is expected

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