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

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90 @ 8637

Last change on this file since 8637 was 8637, checked in by gm, 7 years ago

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8626

File size: 21.0 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         !    fr1_i0  , fr2_i0     = 1sr & 2nd fraction of qsr penetration in ice  [%]
156         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s]
157         !    sprecip              = solid precipitation                           [Kg/m2/s]
158         !    evap_ice             = sublimation                                   [Kg/m2/s]
159         !    qsr_tot , qns_tot    = solar & non solar heat flux (total)           [W/m2]
160         !    qsr_ice , qns_ice    = solar & non solar heat flux over ice          [W/m2]
161         !    dqns_ice             = non solar  heat sensistivity                  [W/m2]
162         !    qemp_oce, qemp_ice,  = sensible heat (associated with evap & precip) [W/m2]
163         !    qprec_ice, qevap_ice
164         !------------------------------------------------------!
165                                        CALL ice_forcing_flx( kt, ksbc )
166         !----------------------------!
167         ! --- ice thermodynamics --- !
168         !----------------------------!
169         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics     
170         !
171         !
172         IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections
173# if defined key_agrif
174         IF( .NOT. Agrif_Root() )       CALL agrif_update_lim3( kt )
175# endif
176         CALL ice_var_glo2eqv        ! necessary calls (at least for coupling)
177         CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling)
178                                     !
179!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
180!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
181!!# if defined key_agrif
182!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid()
183!!# endif
184                                        CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
185!!# if defined key_agrif
186!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid()
187!!# endif
188         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs
189         !
190                                        CALL ice_wri( kt )            ! -- Ice outputs
191         !
192         IF( lrst_ice )                 CALL ice_rst_write( kt )      ! -- Ice restart file
193         !
194         IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash
195         !
196      ENDIF   ! End sea-ice time step only
197
198      !-------------------------!
199      ! --- Ocean time step --- !
200      !-------------------------!
201      IF( ln_icedyn )                   CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses
202!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
203      !
204      IF( nn_timing == 1 )   CALL timing_stop('ice_stp')
205      !
206   END SUBROUTINE ice_stp
207
208
209   SUBROUTINE ice_init
210      !!----------------------------------------------------------------------
211      !!                  ***  ROUTINE ice_init  ***
212      !!
213      !! ** purpose :   Initialize sea-ice parameters
214      !!----------------------------------------------------------------------
215      INTEGER :: ji, jj, ierr
216      !!----------------------------------------------------------------------
217      IF(lwp) WRITE(numout,*)
218      IF(lwp) WRITE(numout,*) 'ice_init: Arrays allocation & Initialization off all routines & init state' 
219      IF(lwp) WRITE(numout,*) '~~~~~~~~'
220      !
221      !                                ! Open the reference and configuration namelist files and namelist output file
222      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
223      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
224      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
225      !
226      CALL par_init                ! set some ice run parameters
227      !
228      !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
229      ierr =        ice_alloc        ()      ! ice variables
230      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
231      ierr = ierr + ice1D_alloc      ()      ! thermodynamics
232      !
233      IF( lk_mpp    )   CALL mpp_sum( ierr )
234      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
235      !
236      CALL ice_itd_init                ! ice thickness distribution initialization
237      !
238      IF( ln_icethd ) THEN
239         CALL ice_thd_init             ! set ice thermodynics parameters (clem: important to call it first for melt ponds)
240      ENDIF   
241      !                                ! Initial sea-ice state
242      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
243         CALL ice_istate_init
244         CALL ice_istate
245      ELSE                                    ! start from a restart file
246         CALL ice_rst_read
247      ENDIF
248      CALL ice_var_glo2eqv
249      CALL ice_var_agg(2)
250      !
251      CALL ice_forcing_init            ! set ice-ocean and ice-atm. coupling parameters
252      !
253      IF( ln_icedyn ) THEN
254         CALL ice_dyn_init             ! set ice dynamics parameters
255      ENDIF
256      !
257      CALL ice_update_init             ! ice surface boundary condition
258      !
259      CALL ice_alb_init                ! ice surface albedo
260      !
261      CALL ice_dia_init                ! initialization for diags
262      !
263      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction
264      tn_ice(:,:,:) = t_su(:,:,:)      ! initialisation of surface temp for coupled simu
265      !
266      !                                ! set max concentration in both hemispheres
267      WHERE( gphit(:,:) > 0._wp )   ;   rn_amax_2d(:,:) = rn_amax_n  ! NH
268      ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH
269      END WHERE
270
271      IF( ln_rstart )   CALL iom_close( numrir )  ! close input ice restart file
272      !
273   END SUBROUTINE ice_init
274
275
276   SUBROUTINE par_init
277      !!-------------------------------------------------------------------
278      !!                  ***  ROUTINE par_init ***
279      !!
280      !! ** Purpose :   Definition generic parameters for ice model
281      !!
282      !! ** Method  :   Read namelist and check the parameter
283      !!                values called at the first timestep (nit000)
284      !!
285      !! ** input   :   Namelist nampar
286      !!-------------------------------------------------------------------
287      INTEGER  ::   ios                 ! Local integer output status for namelist read
288      NAMELIST/nampar/ jpl, nlay_i, nlay_s, nn_monocat, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s,  &
289         &             cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir
290      !!-------------------------------------------------------------------
291      !
292      REWIND( numnam_ice_ref )      ! Namelist nampar in reference namelist : Parameters for ice
293      READ  ( numnam_ice_ref, nampar, IOSTAT = ios, ERR = 901)
294901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampar in reference namelist', lwp )
295
296      REWIND( numnam_ice_cfg )      ! Namelist nampar in configuration namelist : Parameters for ice
297      READ  ( numnam_ice_cfg, nampar, IOSTAT = ios, ERR = 902 )
298902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampar in configuration namelist', lwp )
299      IF(lwm) WRITE ( numoni, nampar )
300      !
301      IF(lwp) THEN                  ! control print
302         WRITE(numout,*)
303         WRITE(numout,*) 'par_init: ice parameters shared among all the routines'
304         WRITE(numout,*) ' ~~~~~~~'
305         WRITE(numout,*) '   Namelist nampar: '
306         WRITE(numout,*) '      number of ice  categories                              jpl    = ', jpl
307         WRITE(numout,*) '      number of ice  layers                                  nlay_i = ', nlay_i
308         WRITE(numout,*) '      number of snow layers                                  nlay_s = ', nlay_s
309         WRITE(numout,*) '      virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat
310         WRITE(numout,*) '      Ice dynamics       (T) or not (F)                   ln_icedyn = ', ln_icedyn
311         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                   ln_icethd = ', ln_icethd
312         WRITE(numout,*) '      maximum ice concentration for NH                              = ', rn_amax_n 
313         WRITE(numout,*) '      maximum ice concentration for SH                              = ', rn_amax_s
314      ENDIF
315      !
316      !                                        !--- check consistency
317      IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN
318         nn_monocat = 0
319         IF(lwp) WRITE(numout,*)
320         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
321      ENDIF
322      IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
323         CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' )
324      ENDIF
325      !
326      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY')
327      !
328      rdt_ice   = REAL(nn_fsbc) * rdt          !--- sea-ice timestep and inverse
329      r1_rdtice = 1._wp / rdt_ice
330      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice = ', rdt_ice
331      !
332      r1_nlay_i = 1._wp / REAL( nlay_i, wp )   !--- inverse of nlay_i and nlay_s
333      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
334      !
335   END SUBROUTINE par_init
336
337
338   SUBROUTINE store_fields
339      !!----------------------------------------------------------------------
340      !!                  ***  ROUTINE store_fields  ***
341      !!
342      !! ** purpose :  store ice variables at "before" time step
343      !!----------------------------------------------------------------------
344      INTEGER  ::   ji, jj, jl      ! dummy loop index
345      !!----------------------------------------------------------------------
346      !
347      a_i_b (:,:,:)   = a_i (:,:,:)     ! ice area
348      v_i_b (:,:,:)   = v_i (:,:,:)     ! ice volume
349      v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume
350      sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content
351      oa_i_b(:,:,:)   = oa_i(:,:,:)     ! areal age content
352      e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy
353      e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy
354      WHERE( a_i_b(:,:,:) >= epsi20 )
355         h_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:)   ! ice thickness
356         h_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:)   ! snw thickness
357      ELSEWHERE
358         h_i_b(:,:,:) = 0._wp
359         h_s_b(:,:,:) = 0._wp
360      END WHERE
361     
362      ! ice velocities & total concentration
363      at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 )
364      u_ice_b(:,:) = u_ice(:,:)
365      v_ice_b(:,:) = v_ice(:,:)
366      !
367   END SUBROUTINE store_fields
368
369
370   SUBROUTINE diag_set0
371      !!----------------------------------------------------------------------
372      !!                  ***  ROUTINE diag_set0  ***
373      !!
374      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
375      !!               of the time step
376      !!----------------------------------------------------------------------
377      INTEGER  ::   ji, jj      ! dummy loop index
378      !!----------------------------------------------------------------------
379      sfx    (:,:) = 0._wp   ;
380      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp
381      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
382      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
383      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
384      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
385      !
386      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
387      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
388      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
389      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
390      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
391      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp 
392      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
393      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
394      wfx_snw_sni(:,:) = 0._wp 
395      wfx_pnd(:,:) = 0._wp
396
397      hfx_thd(:,:) = 0._wp   ;
398      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
399      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
400      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
401      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
402      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
403      hfx_err_rem(:,:) = 0._wp
404      hfx_err_dif(:,:) = 0._wp
405      wfx_err_sub(:,:) = 0._wp
406      !
407      afx_tot(:,:) = 0._wp   ;
408      !
409      diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp
410      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
411
412      ! SIMIP diagnostics
413      diag_fc_bo(:,:)    = 0._wp ; diag_fc_su(:,:)    = 0._wp
414
415      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
416     
417   END SUBROUTINE diag_set0
418
419#else
420   !!----------------------------------------------------------------------
421   !!   Default option           Dummy module         NO ESIM sea-ice model
422   !!----------------------------------------------------------------------
423CONTAINS
424   SUBROUTINE ice_stp ( kt, ksbc )     ! Dummy routine
425      INTEGER, INTENT(in) ::   kt, ksbc
426      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
427   END SUBROUTINE ice_stp
428   SUBROUTINE ice_init                 ! Dummy routine
429   END SUBROUTINE ice_init
430#endif
431
432   !!======================================================================
433END MODULE icestp
Note: See TracBrowser for help on using the repository browser.