source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90 @ 8554

Last change on this file since 8554 was 8534, checked in by clem, 3 years ago

changes in style - part6 - pure cosmetics

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