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

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

Last change on this file was 8892, checked in by frrh, 6 years ago

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