source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90 @ 12475

Last change on this file since 12475 was 12475, checked in by dancopsey, 12 months ago
  • Add more print statements.
  • Move away from using snow to ice diagnostics and use a new snow to pond one instead.
File size: 27.2 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 icecor         ! sea-ice: corrections
58   USE iceupdate      ! sea-ice: sea surface boundary condition update
59   USE icedia         ! sea-ice: budget diagnostics
60   USE icewri         ! sea-ice: outputs
61   USE icerst         ! sea-ice: restarts
62   USE icevar         ! sea-ice: operations
63   USE icectl         ! sea-ice: control
64   USE iceistate      ! sea-ice: initial state
65   USE iceitd         ! sea-ice: remapping thickness distribution
66   USE icealb         ! sea-ice: albedo
67   !
68   USE bdy_oce , ONLY : ln_bdy   ! flag for bdy
69   USE bdyice         ! unstructured open boundary data for sea-ice
70# if defined key_agrif
71   USE agrif_ice
72   USE agrif_ice_interp
73# endif
74   !
75   USE in_out_manager ! I/O manager
76   USE iom            ! I/O manager library
77   USE lib_mpp        ! MPP library
78   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
79   USE timing         ! Timing
80   USE prtctl         ! Print control
81
82   USE lbclnk         ! lateral boundary conditions (or mpp links)
83
84   IMPLICIT NONE
85   PRIVATE
86
87   PUBLIC   ice_stp    ! called by sbcmod.F90
88   PUBLIC   ice_init   ! called by sbcmod.F90
89
90   !! * Substitutions
91#  include "vectopt_loop_substitute.h90"
92   !!----------------------------------------------------------------------
93   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
94   !! $Id$
95   !! Software governed by the CeCILL license (see ./LICENSE)
96   !!----------------------------------------------------------------------
97CONTAINS
98
99   SUBROUTINE ice_stp( kt, ksbc )
100      !!---------------------------------------------------------------------
101      !!                  ***  ROUTINE ice_stp  ***
102      !!
103      !! ** Purpose :   sea-ice model time-stepping and update ocean surface
104      !!              boundary condition over ice-covered area
105      !!
106      !! ** Method  :   ice model time stepping
107      !!              - call the ice dynamics routine
108      !!              - call the ice advection/diffusion routine
109      !!              - call the ice thermodynamics routine
110      !!              - call the routine that computes mass and
111      !!                heat fluxes at the ice/ocean interface
112      !!              - save the outputs
113      !!              - save the outputs for restart when necessary
114      !!
115      !! ** Action  : - time evolution of the LIM sea-ice model
116      !!              - update all sbc variables below sea-ice:
117      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx
118      !!---------------------------------------------------------------------
119      INTEGER, INTENT(in) ::   kt      ! ocean time step
120      INTEGER, INTENT(in) ::   ksbc    ! flux formulation (user defined, bulk, or Pure Coupled)
121      !
122      INTEGER  :: ji, jj, jk, jl, jn, gi, gj   ! dummy loop indices
123      REAL(wp) :: delta_lon, delta_lat, lon_west, lon_east, lat_south, lat_north, lon_print, lat_print
124      !!----------------------------------------------------------------------
125      !
126      IF( ln_timing )   CALL timing_start('ice_stp')
127
128      DO jj = 1, jpj
129         DO ji = 1, jpi
130
131            gi = nimpp + ji - 1
132            gj = njmpp + jj - 1
133            delta_lon = (glamt(2,jj) - glamt(1,jj))/2.0
134            delta_lat = (gphit(ji,2) - gphit(ji,1))/2.0
135
136            lon_west = glamt(ji,jj) - delta_lon
137            lon_east = glamt(ji,jj) + delta_lon
138            lat_south = gphit(ji,jj) - delta_lat
139            lat_north = gphit(ji,jj) + delta_lat
140            to_print_2d(ji,jj) = 0
141
142            lon_print = 95.0
143            lat_print = -63.0  ! For bottom melt investigation
144
145            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
146              to_print_2d(ji,jj) = 1
147              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
148            ENDIF
149
150            lon_print = 10.0
151            lat_print = -62.0   ! For top melt investigation
152
153            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
154              to_print_2d(ji,jj) = 2
155              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
156            ENDIF
157
158            lon_print = 137.81
159            lat_print = -66   ! For conductivity investigation
160
161            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
162              to_print_2d(ji,jj) = 3
163              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
164            ENDIF
165
166            lon_print = 137.81 - 360.0
167
168            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
169              to_print_2d(ji,jj) = 3
170              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
171            ENDIF
172
173            lon_print = 8
174            lat_print = -67   ! For top melt investigation
175
176            IF ( lon_west <= lon_print .AND. lon_print <= lon_east .AND. lat_south <= lat_print .AND. lat_print <= lat_north ) THEN
177              to_print_2d(ji,jj) = 10
178              write(numout,*)'icethd: location to print: ji, jj, gi,gj, to_print_2d(ji,jj) = ', ji, jj, gi,gj, to_print_2d(ji,jj)
179            ENDIF
180          END DO
181        END DO
182
183      !
184      !                                      !-----------------------!
185      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN   ! --- Ice time step --- !
186         !                                   !-----------------------!
187         !
188         kt_ice = kt                              ! -- Ice model time step
189
190         write(numout,*)'icestp 0: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
191         !
192         u_oce(:,:) = ssu_m(:,:)                  ! -- mean surface ocean current
193         v_oce(:,:) = ssv_m(:,:)
194         !
195         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )   ! -- freezing temperature [Kelvin] (set to rt0 over land)
196         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
197         !
198         !                          !==  AGRIF Parent to Child  ==!
199#if defined key_agrif
200         !                              ! nbstep_ice ranges from 1 to the nb of child ocean steps inside one parent ice step
201         IF( .NOT. Agrif_Root() )       nbstep_ice = MOD( nbstep_ice, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
202         !                              ! these calls must remain here for restartability purposes
203                                        CALL agrif_interp_ice( 'T' ) 
204                                        CALL agrif_interp_ice( 'U' )
205                                        CALL agrif_interp_ice( 'V' )
206#endif
207                                        CALL store_fields             ! Store now ice values
208         !
209         !------------------------------------------------!
210         ! --- Dynamical coupling with the atmosphere --- !
211         !------------------------------------------------!
212         ! It provides the following fields used in sea ice model:
213         !    utau_ice, vtau_ice = surface ice stress [N/m2]
214         !------------------------------------------------!
215                                        CALL ice_sbc_tau( kt, ksbc, utau_ice, vtau_ice )     
216
217         write(numout,*)'icestp 1: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
218   
219         !-------------------------------------!
220         ! --- ice dynamics and advection  --- !
221         !-------------------------------------!
222                                        CALL diag_set0                ! set diag of mass, heat and salt fluxes to 0
223                                        CALL ice_rst_opn( kt )        ! Open Ice restart file (if necessary)
224         !
225         IF( ln_icedyn .AND. .NOT.lk_c1d )   &
226            &                           CALL ice_dyn( kt )            ! -- Ice dynamics
227
228         write(numout,*)'icestp 2: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
229         !
230         !                          !==  lateral boundary conditions  ==!
231         IF( ln_icethd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
232
233         write(numout,*)'icestp 3: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
234         !
235         !                          !==  previous lead fraction and ice volume for flux calculations
236                                        CALL ice_var_glo2eqv          ! h_i and h_s for ice albedo calculation
237
238         write(numout,*)'icestp 4: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
239
240                                        CALL ice_var_agg(1)           ! at_i for coupling
241
242         write(numout,*)'icestp 5: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
243
244                                        CALL store_fields             ! Store now ice values
245
246         write(numout,*)'icestp 6: sv_i, e_i, u_ice, e2u, u_oce = ',sv_i(3,4,1), ' ', e_i (3,4,1,1), ' ', u_ice(3,4), ' ', e2u(3,4), ' ', u_oce(3,4)
247         !
248         !------------------------------------------------------!
249         ! --- Thermodynamical coupling with the atmosphere --- !
250         !------------------------------------------------------!
251         ! It provides the following fields used in sea ice model:
252         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s]
253         !    sprecip              = solid precipitation                           [Kg/m2/s]
254         !    evap_ice             = sublimation                                   [Kg/m2/s]
255         !    qsr_tot , qns_tot    = solar & non solar heat flux (total)           [W/m2]
256         !    qsr_ice , qns_ice    = solar & non solar heat flux over ice          [W/m2]
257         !    dqns_ice             = non solar  heat sensistivity                  [W/m2]
258         !    qemp_oce, qemp_ice,  = sensible heat (associated with evap & precip) [W/m2]
259         !    qprec_ice, qevap_ice
260         !------------------------------------------------------!
261                                        CALL ice_sbc_flx( kt, ksbc )
262         write(numout,*)'icestp 7: sv_i, e_s = ',sv_i(3,4,1), ' ', e_s (3,4,1,1)
263         !----------------------------!
264         ! --- ice thermodynamics --- !
265         !----------------------------!
266         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics   
267         write(numout,*)'icestp 8: sv_i, e_s = ',sv_i(3,4,1), ' ', e_s (3,4,1,1)   
268         !
269                                        CALL ice_cor( kt , 2 )        ! -- Corrections
270         write(numout,*)'icestp 9: sv_i, e_s = ',sv_i(3,4,1), ' ', e_s (3,4,1,1)
271         !
272                                        CALL ice_var_glo2eqv          ! necessary calls (at least for coupling)
273         write(numout,*)'icestp 10: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
274                                        CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling)
275         write(numout,*)'icestp 11: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
276         !
277                                        CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
278         write(numout,*)'icestp 12: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
279         !
280         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs
281         !
282                                        CALL ice_wri( kt )            ! -- Ice outputs
283         !
284         IF( lrst_ice )                 CALL ice_rst_write( kt )      ! -- Ice restart file
285         !
286         IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash
287         write(numout,*)'icestp 13: t_s, e_s = ',t_s(ji,jj,jk,jl) - rt0, ' ', e_s (42,26,1,1)
288         !
289      ENDIF   ! End sea-ice time step only
290
291      !-------------------------!
292      ! --- Ocean time step --- !
293      !-------------------------!
294      IF( ln_icedyn )                   CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses
295!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
296      !
297      IF( ln_timing )   CALL timing_stop('ice_stp')
298      !
299   END SUBROUTINE ice_stp
300
301
302   SUBROUTINE ice_init
303      !!----------------------------------------------------------------------
304      !!                  ***  ROUTINE ice_init  ***
305      !!
306      !! ** purpose :   Initialize sea-ice parameters
307      !!----------------------------------------------------------------------
308      INTEGER :: ji, jj, ierr
309      !!----------------------------------------------------------------------
310      IF(lwp) WRITE(numout,*)
311      IF(lwp) WRITE(numout,*) 'Sea Ice Model: SI3 (Sea Ice modelling Integrated Initiative)' 
312      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
313      IF(lwp) WRITE(numout,*)
314      IF(lwp) WRITE(numout,*) 'ice_init: Arrays allocation & Initialization of all routines & init state' 
315      IF(lwp) WRITE(numout,*) '~~~~~~~~'
316      !
317      !                                ! Open the reference and configuration namelist files and namelist output file
318      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
319      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
320      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
321      !
322      CALL par_init                ! set some ice run parameters
323      !
324      !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
325      ierr =        ice_alloc        ()      ! ice variables
326      ierr = ierr + sbc_ice_alloc    ()      ! surface boundary conditions
327      ierr = ierr + ice1D_alloc      ()      ! thermodynamics
328      !
329      CALL mpp_sum( 'icestp', ierr )
330      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
331      !
332      CALL ice_itd_init                ! ice thickness distribution initialization
333      !
334      CALL ice_thd_init                ! set ice thermodynics parameters (clem: important to call it first for melt ponds)
335      !
336      !                                ! Initial sea-ice state
337      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
338         CALL ice_istate_init
339         CALL ice_istate
340      ELSE                                    ! start from a restart file
341         CALL ice_rst_read
342      ENDIF
343      CALL ice_var_glo2eqv
344      CALL ice_var_agg(1)
345      !
346      CALL ice_sbc_init                ! set ice-ocean and ice-atm. coupling parameters
347      !
348      CALL ice_dyn_init                ! set ice dynamics parameters
349      !
350      CALL ice_update_init             ! ice surface boundary condition
351      !
352      CALL ice_alb_init                ! ice surface albedo
353      !
354      CALL ice_dia_init                ! initialization for diags
355      !
356      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction
357      tn_ice(:,:,:) = t_su(:,:,:)      ! initialisation of surface temp for coupled simu
358      !
359      !                                ! set max concentration in both hemispheres
360      WHERE( gphit(:,:) > 0._wp )   ;   rn_amax_2d(:,:) = rn_amax_n  ! NH
361      ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH
362      END WHERE
363
364      IF( ln_rstart )   CALL iom_close( numrir )  ! close input ice restart file
365      !
366   END SUBROUTINE ice_init
367
368
369   SUBROUTINE par_init
370      !!-------------------------------------------------------------------
371      !!                  ***  ROUTINE par_init ***
372      !!
373      !! ** Purpose :   Definition generic parameters for ice model
374      !!
375      !! ** Method  :   Read namelist and check the parameter
376      !!                values called at the first timestep (nit000)
377      !!
378      !! ** input   :   Namelist nampar
379      !!-------------------------------------------------------------------
380      INTEGER  ::   ios                 ! Local integer
381      !!
382      NAMELIST/nampar/ jpl, nlay_i, nlay_s, ln_virtual_itd, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s,  &
383         &             cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir
384      !!-------------------------------------------------------------------
385      !
386      REWIND( numnam_ice_ref )      ! Namelist nampar in reference namelist : Parameters for ice
387      READ  ( numnam_ice_ref, nampar, IOSTAT = ios, ERR = 901)
388901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampar in reference namelist', lwp )
389      REWIND( numnam_ice_cfg )      ! Namelist nampar in configuration namelist : Parameters for ice
390      READ  ( numnam_ice_cfg, nampar, IOSTAT = ios, ERR = 902 )
391902   IF( ios > 0 )   CALL ctl_nam ( ios , 'nampar in configuration namelist', lwp )
392      IF(lwm) WRITE( numoni, nampar )
393      !
394      IF(lwp) THEN                  ! control print
395         WRITE(numout,*)
396         WRITE(numout,*) '   par_init: ice parameters shared among all the routines'
397         WRITE(numout,*) '   ~~~~~~~~'
398         WRITE(numout,*) '      Namelist nampar: '
399         WRITE(numout,*) '         number of ice  categories                           jpl       = ', jpl
400         WRITE(numout,*) '         number of ice  layers                               nlay_i    = ', nlay_i
401         WRITE(numout,*) '         number of snow layers                               nlay_s    = ', nlay_s
402         WRITE(numout,*) '         virtual ITD param for jpl=1 (T) or not (F)     ln_virtual_itd = ', ln_virtual_itd
403         WRITE(numout,*) '         Ice dynamics       (T) or not (F)                   ln_icedyn = ', ln_icedyn
404         WRITE(numout,*) '         Ice thermodynamics (T) or not (F)                   ln_icethd = ', ln_icethd
405         WRITE(numout,*) '         maximum ice concentration for NH                              = ', rn_amax_n 
406         WRITE(numout,*) '         maximum ice concentration for SH                              = ', rn_amax_s
407      ENDIF
408      !                                        !--- change max ice concentration for roundoff errors
409      rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 )
410      rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 )
411      !                                        !--- check consistency
412      IF ( jpl > 1 .AND. ln_virtual_itd ) THEN
413         ln_virtual_itd = .FALSE.
414         IF(lwp) WRITE(numout,*)
415         IF(lwp) WRITE(numout,*) '   ln_virtual_itd forced to false as jpl>1, no need with multiple categories to emulate them'
416      ENDIF
417      !
418      IF( ln_cpl .AND. nn_cats_cpl /= 1 .AND. nn_cats_cpl /= jpl ) THEN
419         CALL ctl_stop( 'STOP', 'par_init: in coupled mode, nn_cats_cpl should be either 1 or jpl' )
420      ENDIF
421      !
422      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY')
423      !
424      rdt_ice   = REAL(nn_fsbc) * rdt          !--- sea-ice timestep and its inverse
425      r1_rdtice = 1._wp / rdt_ice
426      IF(lwp) WRITE(numout,*)
427      IF(lwp) WRITE(numout,*) '      ice timestep rdt_ice = nn_fsbc*rdt = ', rdt_ice
428      !
429      r1_nlay_i = 1._wp / REAL( nlay_i, wp )   !--- inverse of nlay_i and nlay_s
430      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
431      !
432   END SUBROUTINE par_init
433
434
435   SUBROUTINE store_fields
436      !!----------------------------------------------------------------------
437      !!                  ***  ROUTINE store_fields  ***
438      !!
439      !! ** purpose :  store ice variables at "before" time step
440      !!----------------------------------------------------------------------
441      INTEGER  ::   ji, jj, jl      ! dummy loop index
442      !!----------------------------------------------------------------------
443      !
444      a_i_b (:,:,:)   = a_i (:,:,:)     ! ice area
445      v_i_b (:,:,:)   = v_i (:,:,:)     ! ice volume
446      v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume
447      sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content
448      oa_i_b(:,:,:)   = oa_i(:,:,:)     ! areal age content
449      e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy
450      e_i_b (:,:,:,:) = e_i (:,:,:,:)   ! ice thermal energy
451      WHERE( a_i_b(:,:,:) >= epsi20 )
452         h_i_b(:,:,:) = v_i_b(:,:,:) / a_i_b(:,:,:)   ! ice thickness
453         h_s_b(:,:,:) = v_s_b(:,:,:) / a_i_b(:,:,:)   ! snw thickness
454      ELSEWHERE
455         h_i_b(:,:,:) = 0._wp
456         h_s_b(:,:,:) = 0._wp
457      END WHERE
458     
459      WHERE( a_ip(:,:,:) >= epsi20 )
460         h_ip_b(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)   ! ice pond thickness
461      ELSEWHERE
462         h_ip_b(:,:,:) = 0._wp
463      END WHERE
464      !
465      ! ice velocities & total concentration
466      at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 )
467      u_ice_b(:,:) = u_ice(:,:)
468      v_ice_b(:,:) = v_ice(:,:)
469      !
470   END SUBROUTINE store_fields
471
472
473   SUBROUTINE diag_set0
474      !!----------------------------------------------------------------------
475      !!                  ***  ROUTINE diag_set0  ***
476      !!
477      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
478      !!               of the time step
479      !!----------------------------------------------------------------------
480      INTEGER  ::   ji, jj      ! dummy loop index
481      !!----------------------------------------------------------------------
482      sfx    (:,:) = 0._wp   ;
483      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp
484      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
485      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
486      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
487      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
488      sfx_pnd(:,:) = 0._wp
489      !
490      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
491      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
492      wfx_snp(:,:) = 0._wp
493      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
494      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
495      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
496      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp 
497      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
498      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
499      wfx_snw_sni(:,:) = 0._wp 
500      wfx_pnd(:,:) = 0._wp
501
502      hfx_thd(:,:) = 0._wp   ;
503      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
504      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
505      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
506      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
507      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
508      hfx_err_rem(:,:) = 0._wp
509      hfx_err_dif(:,:) = 0._wp
510      wfx_err_sub(:,:) = 0._wp
511      !
512      afx_tot(:,:) = 0._wp   ;
513      !
514      diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp
515      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
516
517      ! SIMIP diagnostics
518      qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes
519      t_si       (:,:,:) = rt0   ! temp at the ice-snow interface
520
521      tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here)
522      cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T)
523      qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T)
524      qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs
525      qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F)
526      !
527      ! for control checks (ln_icediachk)
528      diag_trp_vi(:,:) = 0._wp   ;   diag_trp_vs(:,:) = 0._wp
529      diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp
530      diag_trp_sv(:,:) = 0._wp
531
532   END SUBROUTINE diag_set0
533
534#else
535   !!----------------------------------------------------------------------
536   !!   Default option           Dummy module         NO SI3 sea-ice model
537   !!----------------------------------------------------------------------
538CONTAINS
539   SUBROUTINE ice_stp ( kt, ksbc )     ! Dummy routine
540      INTEGER, INTENT(in) ::   kt, ksbc
541      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
542   END SUBROUTINE ice_stp
543   SUBROUTINE ice_init                 ! Dummy routine
544   END SUBROUTINE ice_init
545#endif
546
547   !!======================================================================
548END MODULE icestp
Note: See TracBrowser for help on using the repository browser.