New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icestp.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

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

changes in style - part5 - start changing init routines

File size: 26.2 KB
RevLine 
[8321]1MODULE icestp
2   !!======================================================================
3   !!                       ***  MODULE  icestp  ***
4   !! Surface module :  update the ocean surface boundary condition over ice
5   !!       &           covered area using LIM sea-ice model
6   !! Sea-Ice model  :  LIM-3 Sea ice model time-stepping
7   !!=====================================================================
8   !! History :  2.0  ! 2006-12  (M. Vancoppenolle) Original code
9   !!            3.0  ! 2008-02  (C. Talandier)  Surface module from icestp.F90
[8411]10   !!             -   ! 2008-04  (G. Madec)  sltyle and ice_ctl routine
[8321]11   !!            3.3  ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
12   !!            3.4  ! 2011-01  (A Porter)  dynamical allocation
[8426]13   !!             -   ! 2012-10  (C. Rousset)  add ice_dia
[8321]14   !!            3.6  ! 2014-07  (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface
15   !!            4.0  ! 2016-06  (L. Brodeau) new unified bulk routine (based on AeroBulk)
16   !!----------------------------------------------------------------------
17#if defined key_lim3
18   !!----------------------------------------------------------------------
[8486]19   !!   'key_lim3'                                    LIM 3.0 sea-ice model
[8321]20   !!----------------------------------------------------------------------
[8486]21   !!   ice_stp       : sea-ice model time-stepping and update ocean surf. boundary cond. over ice-covered area
22   !!   ice_init      :
23   !!   ice_run_init  :
[8321]24   !!----------------------------------------------------------------------
[8486]25   USE oce            ! ocean dynamics and tracers
26   USE dom_oce        ! ocean space and time domain
27   USE c1d            ! 1D vertical configuration
28   USE ice            ! sea-ice: variables
29   USE ice1D          ! sea-ice: thermodynamical 1D variables
[8321]30   !
[8486]31   USE sbc_oce        ! Surface boundary condition: ocean fields
32   USE sbc_ice        ! Surface boundary condition: ice   fields
33   USE iceforcing     ! sea-ice: Surface boundary condition       !!gm why not icesbc module name
[8321]34   !
[8486]35   USE phycst         ! Define parameters for the routines
36   USE eosbn2         ! equation of state
37   USE icerhg         ! sea-ice: rheology
38   USE iceadv         ! sea-ice: advection
39   USE icethd         ! sea-ice: thermodynamics
40   USE icerdgrft      ! sea-ice: ridging/rafting
41   USE iceupdate      ! sea-ice: sea surface boundary condition update
42   USE icedia         ! sea-ice: budget diagnostics
43   USE icewri         ! sea-ice: outputs
44   USE icerst         ! sea-ice: restarts
45   USE icecor         ! sea-ice: corrections
46   USE icevar         ! sea-ice: operations
47   USE icectl         ! sea-ice: control
[8321]48   ! MV MP 2016
[8486]49   USE limmp          ! sea-ice: melt ponds
[8321]50   ! END MV MP 2016
[8486]51   USE iceist         ! sea-ice: initial state
52   USE icethd_sal     ! sea-ice: thermodynamics and salinity
[8505]53   USE iceitd         ! sea-ice: remapping thickness distribution
54   USE icealb         ! sea-ice: albedo
[8321]55   !
[8486]56   USE bdy_oce , ONLY : ln_bdy   ! flag for bdy
57   USE bdyice         ! unstructured open boundary data for sea-ice
[8321]58# if defined key_agrif
59   USE agrif_ice
60   USE agrif_lim3_update
61   USE agrif_lim3_interp
62# endif
[8486]63   !
64   USE in_out_manager ! I/O manager
65   USE iom            ! I/O manager library
66   USE prtctl         ! Print control
67   USE lib_fortran    !
68   USE lbclnk         ! lateral boundary condition - MPP link
69   USE lib_mpp        ! MPP library
70   USE timing         ! Timing
[8321]71
72   IMPLICIT NONE
73   PRIVATE
74
[8486]75   PUBLIC   ice_stp    ! called by sbcmod.F90
76   PUBLIC   ice_init   ! called by sbcmod.F90
[8321]77
[8486]78   INTEGER ::              nice_dyn   ! choice of the type of advection scheme
79   !                                  ! associated indices:
80   INTEGER, PARAMETER ::   np_dynNO   = 0   ! no ice dynamics and ice advection
81   INTEGER, PARAMETER ::   np_dynFULL = 1   ! full ice dynamics  (rheology + advection + ridging/rafting + correction)
82   INTEGER, PARAMETER ::   np_dyn     = 2   ! no ridging/rafting (rheology + advection                   + correction)
83   INTEGER, PARAMETER ::   np_dynPURE = 3   ! pure dynamics      (rheology + advection)
84
[8321]85   !! * Substitutions
86#  include "vectopt_loop_substitute.h90"
87   !!----------------------------------------------------------------------
[8486]88   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
[8321]89   !! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $
90   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
91   !!----------------------------------------------------------------------
92CONTAINS
93
94   SUBROUTINE ice_stp( kt, ksbc )
95      !!---------------------------------------------------------------------
96      !!                  ***  ROUTINE ice_stp  ***
97      !!
[8486]98      !! ** Purpose :   sea-ice model time-stepping and update ocean surface
99      !!              boundary condition over ice-covered area
[8321]100      !!
101      !! ** Method  :   ice model time stepping
102      !!              - call the ice dynamics routine
103      !!              - call the ice advection/diffusion routine
104      !!              - call the ice thermodynamics routine
105      !!              - call the routine that computes mass and
106      !!                heat fluxes at the ice/ocean interface
107      !!              - save the outputs
108      !!              - save the outputs for restart when necessary
109      !!
110      !! ** Action  : - time evolution of the LIM sea-ice model
111      !!              - update all sbc variables below sea-ice:
112      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx
113      !!---------------------------------------------------------------------
114      INTEGER, INTENT(in) ::   kt      ! ocean time step
[8486]115      INTEGER, INTENT(in) ::   ksbc    ! flux formulation (user defined, bulk, or Pure Coupled)
116      !
117      INTEGER ::   jl   ! dummy loop index
[8321]118      !!----------------------------------------------------------------------
[8486]119      !
[8321]120      IF( nn_timing == 1 )   CALL timing_start('ice_stp')
[8486]121      !
122      !                                      !-----------------------!
123      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN   ! --- Ice time step --- !
124         !                                   !-----------------------!
125         !
[8498]126         kt_ice = kt                              ! -- Ice model time step
[8486]127         !
[8321]128# if defined key_agrif
129         IF( .NOT. Agrif_Root() )  lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
130# endif
[8498]131         u_oce(:,:) = ssu_m(:,:)                  ! -- mean surface ocean current
132         v_oce(:,:) = ssv_m(:,:)
133         !
134         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )   ! -- freezing temperature [Kelvin] (set to rt0 over land)
[8321]135         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
136         !
[8498]137         CALL store_fields                        ! -- Store now ice values
[8486]138         !
[8321]139         !------------------------------------------------!
140         ! --- Dynamical coupling with the atmosphere --- !
141         !------------------------------------------------!
[8486]142         ! It provides the following fields used in sea ice model:
[8404]143         !    utau_ice, vtau_ice = surface ice stress [N/m2]
[8486]144         !------------------------------------------------!
145         CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )
[8404]146                                     
[8486]147         !-------------------------------------!
148         ! --- ice dynamics and advection  --- !
149         !-------------------------------------!
150         CALL ice_diag0             ! set diag of mass, heat and salt fluxes to 0
151         CALL ice_rst_opn( kt )     ! Open Ice restart file (if necessary)
[8321]152         !
[8486]153         SELECT CASE( nice_dyn )   
154         CASE ( np_dynFULL )          !==  all dynamical processes  ==!
155            CALL ice_rhg   ( kt )          ! -- rheology 
156            CALL ice_adv   ( kt )          ! -- advection of ice
157            CALL ice_rdgrft( kt )          ! -- ridging/rafting
158            CALL ice_cor   ( kt , 1 )      ! -- Corrections
159         CASE ( np_dyn )              !==  pure dynamics only ==!   (no ridging/rafting)   (nono cat. case 2)
160            CALL ice_rhg   ( kt )          ! -- rheology 
161            CALL ice_adv   ( kt )          ! -- advection of ice
162            CALL ice_cor   ( kt , 1 )      ! -- Corrections
163         CASE ( np_dynPURE )          !==  pure dynamics only ==!   (nn_limdyn= 0 or 1 )
164            CALL ice_rhg   ( kt )          ! -- rheology 
165            CALL ice_adv   ( kt )          ! -- advection of ice
166         END SELECT
167         !
168         !                          !==  lateral boundary conditions  ==!
[8321]169#if defined key_agrif
[8486]170         IF( .NOT. Agrif_Root()     )   CALL agrif_interp_lim3('T')   ! -- AGRIF
[8321]171#endif
[8486]172         IF( ln_limthd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
173         !
174         !
175         !                          !==  previous lead fraction and ice volume for flux calculations
176         !
177         CALL ice_var_glo2eqv            ! ht_i and ht_s for ice albedo calculation
178         CALL ice_var_agg(1)             ! at_i for coupling
[8498]179         CALL store_fields               ! Store now ice values
[8321]180
181         !------------------------------------------------------!
182         ! --- Thermodynamical coupling with the atmosphere --- !
183         !------------------------------------------------------!
[8404]184         ! It provides the following fields used in sea ice model:
[8486]185         !    fr1_i0  , fr2_i0     = 1sr & 2nd fraction of qsr penetration in ice  [%]
186         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s]
187         !    sprecip              = solid precipitation                           [Kg/m2/s]
188         !    evap_ice             = sublimation                                   [Kg/m2/s]
189         !    qsr_tot , qns_tot    = solar & non solar heat flux (total)           [W/m2]
190         !    qsr_ice , qns_ice    = solar & non solar heat flux over ice          [W/m2]
191         !    dqns_ice             = non solar  heat sensistivity                  [W/m2]
192         !    qemp_oce, qemp_ice,  = sensible heat (associated with evap & precip) [W/m2]
193         !    qprec_ice, qevap_ice
194         !------------------------------------------------------!
195         CALL ice_forcing_flx( kt, ksbc )
[8321]196
197         !----------------------------!
198         ! --- ice thermodynamics --- !
199         !----------------------------!
[8486]200         IF( ln_limthd )            CALL ice_thd( kt )          ! -- Ice thermodynamics     
[8321]201
202         ! MV MP 2016
[8486]203         IF ( ln_pnd )              CALL lim_mp( kt )           ! -- Melt ponds
[8321]204         ! END MV MP 2016
205
[8486]206         IF( ln_limthd )            CALL ice_cor( kt , 2 )      ! -- Corrections
[8321]207         ! ---
208# if defined key_agrif
[8486]209         IF( .NOT. Agrif_Root() )   CALL agrif_update_lim3( kt )
[8321]210# endif
[8486]211                                    CALL ice_var_glo2eqv        ! necessary calls (at least for coupling)
212                                    CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling)
213                                    !
[8321]214!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
215!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
216!!# if defined key_agrif
217!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid()
218!!# endif
[8486]219                                    CALL ice_update_flx( kt )   ! -- Update ocean surface mass, heat and salt fluxes
[8321]220!!# if defined key_agrif
221!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
222!!# endif
[8498]223         IF( ln_limdiahsb )         CALL ice_dia( kt )        ! -- Diagnostics and outputs
[8321]224         !
[8498]225                                    CALL ice_wri( 1 )         ! -- Ice outputs
[8321]226         !
[8486]227         IF( kt == nit000 .AND. ln_rstart )   &                !!gm  I don't understand the ln_rstart, if needed, add a comment, please
[8498]228            &                       CALL iom_close( numrir )  ! close input ice restart file
[8321]229         !
[8498]230         IF( lrst_ice )             CALL ice_rst_write( kt )  ! -- Ice restart file
[8321]231         !
[8498]232         IF( ln_limctl )            CALL ice_ctl( kt )        ! alerts in case of model crash
[8321]233         !
234      ENDIF   ! End sea-ice time step only
235
236      !-------------------------!
237      ! --- Ocean time step --- !
238      !-------------------------!
239      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
240      !    using before instantaneous surf. currents
[8498]241      IF( ln_limdyn )               CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )
[8321]242!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
243      !
244      IF( nn_timing == 1 ) CALL timing_stop('ice_stp')
245      !
246   END SUBROUTINE ice_stp
247
248
249   SUBROUTINE ice_init
250      !!----------------------------------------------------------------------
251      !!                  ***  ROUTINE ice_init  ***
252      !!
253      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules
254      !!----------------------------------------------------------------------
255      INTEGER :: ji, jj, ierr
256      !!----------------------------------------------------------------------
257      IF(lwp) WRITE(numout,*)
258      IF(lwp) WRITE(numout,*) 'ice_init : update ocean surface boudary condition' 
259      IF(lwp) WRITE(numout,*) '~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
260      !
261      !                                ! Open the reference and configuration namelist files and namelist output file
262      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
263      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
264      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
265      !
266      CALL ice_run_init                ! set some ice run parameters
267      !
[8324]268      !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
[8321]269      ierr =        ice_alloc        ()      ! ice variables
[8331]270      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
[8420]271      ierr = ierr + ice1D_alloc      ()      ! thermodynamics
[8409]272      IF( ln_limdyn )   ierr = ierr + ice_rdgrft_alloc ()      ! ridging/rafting
[8321]273      !
274      IF( lk_mpp    )   CALL mpp_sum( ierr )
275      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
276      !
[8407]277      CALL ice_rhg_init                ! set ice dynamics parameters
[8321]278      !
279      CALL ice_itd_init                ! ice thickness distribution initialization
280      !
[8422]281      CALL ice_thd_init                ! set ice thermodynics parameters
[8321]282      !
[8422]283      CALL ice_thd_sal_init            ! set ice salinity parameters
[8321]284       
285      ! MV MP 2016
286      CALL lim_mp_init                 ! set melt ponds parameters
287      ! END MV MP 2016
288
[8505]289      IF( ln_limdyn )   CALL ice_rdgrft_init  ! ice thickness distribution initialization for ridging/rafting
[8321]290      !                                ! Initial sea-ice state
291      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
[8426]292         CALL ice_ist
[8321]293      ELSE                                    ! start from a restart file
[8413]294         CALL ice_rst_read
[8321]295      ENDIF
[8424]296      CALL ice_var_agg(2)
297      CALL ice_var_glo2eqv
[8321]298      !
[8414]299      CALL ice_update_init                 ! ice surface boundary condition
[8321]300      !
[8505]301      CALL ice_alb_init                    ! ice surface albedo
[8321]302      !
[8505]303      CALL ice_dia_init                    ! initialization for diags
304      !
[8486]305      fr_i  (:,:)   = at_i(:,:)         ! initialisation of sea-ice fraction
[8321]306      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
307      !
308      DO jj = 1, jpj
309         DO ji = 1, jpi
[8486]310            IF( gphit(ji,jj) > 0._wp ) THEN   ;   rn_amax_2d(ji,jj) = rn_amax_n  ! NH
311            ELSE                              ;   rn_amax_2d(ji,jj) = rn_amax_s  ! SH
[8321]312            ENDIF
313         END DO
314      END DO
315      !
316   END SUBROUTINE ice_init
317
318
319   SUBROUTINE ice_run_init
320      !!-------------------------------------------------------------------
321      !!                  ***  ROUTINE ice_run_init ***
322      !!
323      !! ** Purpose :   Definition some run parameter for ice model
324      !!
325      !! ** Method  :   Read the namicerun namelist and check the parameter
326      !!              values called at the first timestep (nit000)
327      !!
328      !! ** input   :   Namelist namicerun
329      !!-------------------------------------------------------------------
330      INTEGER  ::   ios                 ! Local integer output status for namelist read
331      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   &
332         &                cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice 
333      !!-------------------------------------------------------------------
334      !
[8486]335      REWIND( numnam_ice_ref )      ! Namelist namicerun in reference namelist : Parameters for ice
[8321]336      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
337901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
338
[8486]339      REWIND( numnam_ice_cfg )      ! Namelist namicerun in configuration namelist : Parameters for ice
[8321]340      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
341902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
342      IF(lwm) WRITE ( numoni, namicerun )
343      !
[8486]344      IF(lwp) THEN                  ! control print
[8321]345         WRITE(numout,*)
346         WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
347         WRITE(numout,*) ' ~~~~~~'
[8486]348         WRITE(numout,*) '   Namelist namicerun : '
349         WRITE(numout,*) '      number of ice  categories                              jpl    = ', jpl
350         WRITE(numout,*) '      number of ice  layers                                  nlay_i = ', nlay_i
351         WRITE(numout,*) '      number of snow layers                                  nlay_s = ', nlay_s
352         WRITE(numout,*) '      virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat
353         WRITE(numout,*) '      maximum ice concentration for NH                              = ', rn_amax_n 
354         WRITE(numout,*) '      maximum ice concentration for SH                              = ', rn_amax_s
355         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                  ln_limthd  = ', ln_limthd
356         WRITE(numout,*) '      Ice dynamics       (T) or not (F)                  ln_limdyn  = ', ln_limdyn
357         WRITE(numout,*) '         associated switch                               nn_limdyn  = ', nn_limdyn
358         WRITE(numout,*) '            =2 all processes (default option)'
359         WRITE(numout,*) '            =1 advection only (no ridging/rafting)'
360         WRITE(numout,*) '            =0 advection only with prescribed velocity given by '
361         WRITE(numout,*) '               a uniform field        (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')'
[8321]362      ENDIF
363      !
[8486]364      !                                         !--- check consistency
365      IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN
[8321]366         nn_monocat = 0
367         IF(lwp) WRITE(numout,*)
368         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
369      ENDIF
[8486]370      IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
[8321]371         CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
372      ENDIF
373      !
[8486]374      IF( ln_bdy .AND. ln_limdiachk )   CALL ctl_warn('ice_run_init: online conservation check does not work with BDY')
[8321]375      !
[8486]376      !                             ! set the choice of ice dynamics
377      IF( lk_c1d .OR. .NOT.ln_limdyn ) THEN
378         nice_dyn = np_dynNO                    !--- no dynamics
379      ELSE
380         SELECT CASE( nn_limdyn )
381         CASE( 2 )                   
382            IF( nn_monocat /= 2 ) THEN          !--- full dynamics (rheology + advection + ridging/rafting + correction)
383               nice_dyn = np_dynFULL
384            ELSE
385               nice_dyn = np_dyn                !--- dynamics without ridging/rafting
386            ENDIF
387         CASE( 0 , 1 )                          !--- dynamics without ridging/rafting and correction
388            nice_dyn = np_dynPURE
389         END SELECT
390      ENDIF
391      !                                         !--- simple conservative piling, comparable with LIM2
392      l_piling = nn_limdyn == 1 .OR. ( nn_monocat == 2  .AND.  jpl == 1 )
[8321]393      !
[8486]394      rdt_ice   = REAL(nn_fsbc) * rdt           !--- sea-ice timestep and inverse
395      r1_rdtice = 1._wp / rdt_ice
[8321]396      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice
397      !
[8486]398      r1_nlay_i = 1._wp / REAL( nlay_i, wp )    !--- inverse of nlay_i and nlay_s
399      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
400      !
[8321]401   END SUBROUTINE ice_run_init
402
403
[8498]404   SUBROUTINE store_fields
[8321]405      !!----------------------------------------------------------------------
[8498]406      !!                  ***  ROUTINE store_fields  ***
[8321]407      !!
408      !! ** purpose :  store ice variables at "before" time step
409      !!----------------------------------------------------------------------
[8360]410      INTEGER  ::   ji, jj, jl      ! dummy loop index
411      !!----------------------------------------------------------------------
[8321]412      !
[8360]413      DO jl = 1, jpl
[8404]414         DO jj = 2, jpjm1
415            DO ji = 2, jpim1
[8360]416               a_i_b  (ji,jj,jl)   = a_i  (ji,jj,jl)     ! ice area
417               v_i_b  (ji,jj,jl)   = v_i  (ji,jj,jl)     ! ice volume
418               v_s_b  (ji,jj,jl)   = v_s  (ji,jj,jl)     ! snow volume
419               smv_i_b(ji,jj,jl)   = smv_i(ji,jj,jl)     ! salt content
420               oa_i_b (ji,jj,jl)   = oa_i (ji,jj,jl)     ! areal age content
421               e_s_b  (ji,jj,:,jl) = e_s  (ji,jj,:,jl)   ! snow thermal energy
422               e_i_b  (ji,jj,:,jl) = e_i  (ji,jj,:,jl)   ! ice thermal energy
423               !                                         ! ice thickness
424               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes
425               ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
426               ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
427            END DO
[8404]428         END DO   
[8360]429      END DO
[8486]430      CALL lbc_lnk_multi(  a_i_b, 'T', 1., v_i_b , 'T', 1., ht_i_b, 'T', 1., smv_i_b, 'T', 1.,   &
431         &                oa_i_b, 'T', 1., v_s_b , 'T', 1., ht_s_b, 'T', 1. )
[8404]432      CALL lbc_lnk( e_i_b, 'T', 1. )
433      CALL lbc_lnk( e_s_b, 'T', 1. )
[8321]434     
[8486]435!!gm Question:  here , a_i_b, u_ice and v_ice  are defined over the whole domain,
436!!gm            so why not just a copy over the whole domain and no lbc_lnk ?
437!!gm            that is some thing like:
438!            at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 )
439!            u_ice_b(:,:) = u_ice(:,:)
440!            v_ice_b(:,:) = v_ice(:,:)
441!    idem for the loop above
442!!gm
[8360]443      ! ice velocities & total concentration
[8404]444      DO jj = 2, jpjm1
445         DO ji = 2, jpim1
[8360]446            at_i_b(ji,jj)  = SUM( a_i_b(ji,jj,:) )
447            u_ice_b(ji,jj) = u_ice(ji,jj)
448            v_ice_b(ji,jj) = v_ice(ji,jj)
449         END DO
450      END DO
[8404]451      CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. )
[8486]452      !
[8498]453   END SUBROUTINE store_fields
[8321]454
455
456   SUBROUTINE ice_diag0
457      !!----------------------------------------------------------------------
458      !!                  ***  ROUTINE ice_diag0  ***
459      !!
460      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
461      !!               of the time step
462      !!----------------------------------------------------------------------
[8360]463      INTEGER  ::   ji, jj      ! dummy loop index
464      !!----------------------------------------------------------------------
465      DO jj = 1, jpj
466         DO ji = 1, jpi
467            sfx    (ji,jj) = 0._wp   ;
468            sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp
469            sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp
470            sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp
471            sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp
472            sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp
473            !
474            wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp
475            wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp
476            wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp
477            wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp
478            wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp
479            wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp 
480            wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
481            wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
482            wfx_snw_sni(ji,jj) = 0._wp 
483            ! MV MP 2016
484            wfx_pnd(ji,jj) = 0._wp
485            ! END MV MP 2016
486           
487            hfx_thd(ji,jj) = 0._wp   ;
488            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
489            hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp
490            hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp
491            hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp
492            hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp
493            hfx_err(ji,jj) = 0._wp   ;   hfx_err_rem(ji,jj) = 0._wp
494            hfx_err_dif(ji,jj) = 0._wp
495            wfx_err_sub(ji,jj) = 0._wp
496            !
497            afx_tot(ji,jj) = 0._wp   ;
498            !
499            diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp
500            diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp
501           
502            ! SIMIP diagnostics
503            diag_fc_bo(ji,jj)    = 0._wp ; diag_fc_su(ji,jj)    = 0._wp
504           
505            tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
506         END DO
507      END DO
[8321]508     
509   END SUBROUTINE ice_diag0
510
511#else
512   !!----------------------------------------------------------------------
513   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
514   !!----------------------------------------------------------------------
515CONTAINS
516   SUBROUTINE ice_stp ( kt, ksbc )     ! Dummy routine
517      INTEGER, INTENT(in) ::   kt, ksbc
518      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
519   END SUBROUTINE ice_stp
520   SUBROUTINE ice_init                 ! Dummy routine
521   END SUBROUTINE ice_init
522#endif
523
524   !!======================================================================
525END MODULE icestp
Note: See TracBrowser for help on using the repository browser.