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 NEMO/branches/UKMO/NEMO_4.0.4_remade_heat_balance/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_remade_heat_balance/src/ICE/icestp.F90 @ 15814

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