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

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

first step to make melt ponds compliant with the new code

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