source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 8291

Last change on this file since 8291 was 8291, checked in by clem, 4 years ago

correct last commit changes. This is the start version for the modifications of the ice model. All sette tests passed except restartability for agrif and sas (though sas_biper is restartable). Configurations CREG, SPITZ, FER & SASBIPER are running.

  • Property svn:keywords set to Id
File size: 34.9 KB
Line 
1MODULE sbcice_lim
2   !!======================================================================
3   !!                       ***  MODULE  sbcice_lim  ***
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
10   !!             -   ! 2008-04  (G. Madec)  sltyle and lim_ctl routine
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
13   !!             -   ! 2012-10  (C. Rousset)  add lim_diahsb
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   !!----------------------------------------------------------------------
19   !!   'key_lim3' :                                  LIM 3.0 sea-ice model
20   !!----------------------------------------------------------------------
21   !!   sbc_ice_lim  : sea-ice model time-stepping and update ocean sbc over ice-covered area
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers
24   USE dom_oce         ! ocean space and time domain
25   USE ice             ! LIM-3: ice variables
26   USE thd_ice         ! LIM-3: thermodynamical variables
27   !
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE sbc_ice         ! Surface boundary condition: ice   fields
30   USE usrdef_sbc      ! user defined: surface boundary condition
31   USE sbcblk          ! Surface boundary condition: bulk
32   USE sbccpl          ! Surface boundary condition: coupled interface
33   USE albedo          ! ocean & ice albedo
34   !
35   USE phycst          ! Define parameters for the routines
36   USE eosbn2          ! equation of state
37   USE limdyn          ! Ice dynamics
38   USE limtrp          ! Ice transport
39   USE limhdf          ! Ice horizontal diffusion
40   USE limthd          ! Ice thermodynamics
41   USE limitd_me       ! Mechanics on ice thickness distribution
42   USE limsbc          ! sea surface boundary condition
43   USE limdiahsb       ! Ice budget diagnostics
44   USE limwri          ! Ice outputs
45   USE limrst          ! Ice restarts
46   USE limupdate1      ! update of global variables
47   USE limupdate2      ! update of global variables
48   USE limvar          ! Ice variables switch
49   USE limctl          !
50   ! MV MP 2016
51   USE limmp
52   ! END MV MP 2016
53
54   USE limistate       ! LIM initial state
55   USE limthd_sal      ! LIM ice thermodynamics: salinity
56   !
57   USE c1d             ! 1D vertical configuration
58   USE in_out_manager  ! I/O manager
59   USE iom             ! I/O manager library
60   USE prtctl          ! Print control
61   USE lib_fortran     !
62   USE lbclnk          ! lateral boundary condition - MPP link
63   USE lib_mpp         ! MPP library
64   USE wrk_nemo        ! work arrays
65   USE timing          ! Timing
66
67   USE bdy_oce   , ONLY: ln_bdy
68   USE bdyice_lim      ! unstructured open boundary data  (bdy_ice_lim routine)
69# if defined key_agrif
70   USE agrif_ice
71   USE agrif_lim3_update
72   USE agrif_lim3_interp
73# endif
74
75   IMPLICIT NONE
76   PRIVATE
77
78   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90
79   PUBLIC sbc_lim_init ! routine called by sbcmod.F90
80
81   !! * Substitutions
82#  include "vectopt_loop_substitute.h90"
83   !!----------------------------------------------------------------------
84   !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
85   !! $Id$
86   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
87   !!----------------------------------------------------------------------
88CONTAINS
89
90   SUBROUTINE sbc_ice_lim( kt, ksbc )
91      !!---------------------------------------------------------------------
92      !!                  ***  ROUTINE sbc_ice_lim  ***
93      !!
94      !! ** Purpose :   update the ocean surface boundary condition via the
95      !!                Louvain la Neuve Sea Ice Model time stepping
96      !!
97      !! ** Method  :   ice model time stepping
98      !!              - call the ice dynamics routine
99      !!              - call the ice advection/diffusion routine
100      !!              - call the ice thermodynamics routine
101      !!              - call the routine that computes mass and
102      !!                heat fluxes at the ice/ocean interface
103      !!              - save the outputs
104      !!              - save the outputs for restart when necessary
105      !!
106      !! ** Action  : - time evolution of the LIM sea-ice model
107      !!              - update all sbc variables below sea-ice:
108      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx
109      !!---------------------------------------------------------------------
110      INTEGER, INTENT(in) ::   kt      ! ocean time step
111      INTEGER, INTENT(in) ::   ksbc    ! type of sbc flux ( 1 = user defined formulation,
112                                       !                    3 = bulk formulation,
113                                       !                    4 = Pure Coupled formulation)
114      !!
115      INTEGER  ::   jl                 ! dummy loop index
116      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
117      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice 
118      !!----------------------------------------------------------------------
119
120      IF( nn_timing == 1 )   CALL timing_start('sbc_ice_lim')
121
122      !-----------------------!
123      ! --- Ice time step --- !
124      !-----------------------!
125      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
126
127# if defined key_agrif
128         IF( .NOT. Agrif_Root() )  lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
129# endif
130
131         ! mean surface ocean current at ice velocity point (C-grid dynamics :  U- & V-points as the ocean)
132         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
133         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
134
135         ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
136         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
137         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
138
139         ! Mask sea ice surface temperature (set to rt0 over land)
140         DO jl = 1, jpl
141            t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
142         END DO
143         !
144         !------------------------------------------------!
145         ! --- Dynamical coupling with the atmosphere --- !
146         !------------------------------------------------!
147         ! It provides the following fields:
148         ! utau_ice, vtau_ice : surface ice stress (U- & V-points)   [N/m2]
149         !-----------------------------------------------------------------
150                                      CALL sbc_lim_bef         ! Store previous ice values
151         SELECT CASE( ksbc )
152            CASE( jp_usr     )   ;    CALL usrdef_sbc_ice_tau( kt )                 ! user defined formulation
153            CASE( jp_blk     )   ;    CALL blk_ice_tau                              ! Bulk formulation
154            CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled   formulation
155         END SELECT
156
157         IF( ln_mixcpl) THEN                                                       ! Case of a mixed Bulk/Coupled formulation
158            CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice)
159                                      CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
160            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
161            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
162            CALL wrk_dealloc( jpi,jpj  , zutau_ice, zvtau_ice)
163         ENDIF
164
165         !-------------------------------------------------------!
166         ! --- ice dynamics and transport (except in 1D case) ---!
167         !-------------------------------------------------------!
168         numit = numit + nn_fsbc                  ! Ice model time step
169         !
170                                      CALL sbc_lim_diag0       ! set diag of mass, heat and salt fluxes to 0
171                                      CALL lim_rst_opn( kt )   ! Open Ice restart file
172         !
173         ! --- zap this if no ice dynamics --- !
174         IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN
175            !
176            IF( nn_limdyn /= 0 ) THEN                          ! -- Ice dynamics
177                                      CALL lim_dyn( kt )       !     rheology 
178            ELSE
179               u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity
180               v_ice(:,:) = rn_vice * vmask(:,:,1)
181               !!CALL RANDOM_NUMBER(u_ice(:,:))
182               !!CALL RANDOM_NUMBER(v_ice(:,:))
183            ENDIF
184                                      CALL lim_trp( kt )       ! -- Ice transport (Advection/diffusion)
185            IF( nn_limdyn == 2 .AND. nn_monocat /= 2 )  &      ! -- Mechanical redistribution (ridging/rafting)
186               &                      CALL lim_itd_me         
187            IF( nn_limdyn == 2 )      CALL lim_update1( kt )   ! -- Corrections
188            !
189         ENDIF
190
191         ! ---
192#if defined key_agrif
193         IF( .NOT. Agrif_Root() )     CALL agrif_interp_lim3('T')
194#endif
195         IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice_lim( kt )   ! -- bdy ice thermo
196         ! previous lead fraction and ice volume for flux calculations
197                                      CALL sbc_lim_bef                       
198                                      CALL lim_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation
199                                      CALL lim_var_agg(1)      ! at_i for coupling (via pfrld)
200         !
201         pfrld(:,:)   = 1._wp - at_i(:,:)
202         phicif(:,:)  = vt_i(:,:)
203
204         !------------------------------------------------------!
205         ! --- Thermodynamical coupling with the atmosphere --- !
206         !------------------------------------------------------!
207         ! It provides the following fields:
208         ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2]
209         ! qla_ice            : latent heat flux over ice              (T-point)         [W/m2]
210         ! dqns_ice, dqla_ice : non solar & latent heat sensistivity   (T-point)         [W/m2]
211         ! tprecip , sprecip  : total & solid precipitation            (T-point)         [Kg/m2/s]
212         ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%]
213         !----------------------------------------------------------------------------------------
214         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs )
215         
216                                      CALL albedo_ice( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos MV MP 2016
217
218         SELECT CASE( ksbc )
219            CASE( jp_usr )   ;        CALL usrdef_sbc_ice_flx( kt ) ! user defined formulation
220            CASE( jp_blk )                                          ! bulk formulation
221               ! albedo depends on cloud fraction because of non-linear spectral effects
222               alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
223                                      CALL blk_ice_flx( t_su, alb_ice )
224               IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
225               IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
226            CASE ( jp_purecpl )
227               ! albedo depends on cloud fraction because of non-linear spectral effects
228               alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
229                                      CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
230               IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
231         END SELECT
232
233         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs )
234
235         !----------------------------!
236         ! --- ice thermodynamics --- !
237         !----------------------------!
238         ! --- zap this if no ice thermo --- !
239         IF( ln_limthd )              CALL lim_thd( kt )        ! -- Ice thermodynamics     
240
241         ! MV MP 2016
242         IF ( ln_pnd )                CALL lim_mp( kt )         ! -- Melt ponds
243         ! END MV MP 2016
244
245         IF( ln_limthd )              CALL lim_update2( kt )    ! -- Corrections
246         ! ---
247# if defined key_agrif
248         IF( .NOT. Agrif_Root() )     CALL agrif_update_lim3( kt )
249# endif
250                                      CALL lim_var_glo2eqv      ! necessary calls (at least for coupling)
251                                      CALL lim_var_agg( 2 )     ! necessary calls (at least for coupling)
252                                      !
253!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
254!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
255!!# if defined key_agrif
256!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid()
257!!# endif
258                                      CALL lim_sbc_flx( kt )    ! -- Update surface ocean mass, heat and salt fluxes
259!!# if defined key_agrif
260!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
261!!# endif
262         IF( ln_limdiahsb )           CALL lim_diahsb( kt )     ! -- Diagnostics and outputs
263         !
264                                      CALL lim_wri( 1 )         ! -- Ice outputs
265         !
266         IF( kt == nit000 .AND. ln_rstart )   &
267            &                         CALL iom_close( numrir )  ! close input ice restart file
268         !
269         IF( lrst_ice )               CALL lim_rst_write( kt )  ! -- Ice restart file
270         !
271         IF( ln_limctl )              CALL lim_ctl( kt )        ! alerts in case of model crash
272         !
273      ENDIF   ! End sea-ice time step only
274
275      !-------------------------!
276      ! --- Ocean time step --- !
277      !-------------------------!
278      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
279      !    using before instantaneous surf. currents
280      IF( ln_limdyn )                 CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )
281!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
282      !
283      IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim')
284      !
285   END SUBROUTINE sbc_ice_lim
286
287
288   SUBROUTINE sbc_lim_init
289      !!----------------------------------------------------------------------
290      !!                  ***  ROUTINE sbc_lim_init  ***
291      !!
292      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules
293      !!----------------------------------------------------------------------
294      INTEGER :: ji, jj, ierr
295      !!----------------------------------------------------------------------
296      IF(lwp) WRITE(numout,*)
297      IF(lwp) WRITE(numout,*) 'sbc_lim_init : update ocean surface boudary condition' 
298      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
299      !
300      !                                ! Open the reference and configuration namelist files and namelist output file
301      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
302      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
303      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
304      !
305      CALL lim_run_init                ! set some ice run parameters
306      !
307      !                                ! Allocate the ice arrays
308      ierr =        ice_alloc        ()      ! ice variables
309      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing
310      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics
311      IF( ln_limdyn )   ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics
312      !
313      IF( lk_mpp    )   CALL mpp_sum( ierr )
314      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'sbc_lim_init : unable to allocate ice arrays')
315      !
316      CALL lim_dyn_init                ! set ice dynamics parameters
317      !
318      CALL lim_itd_init                ! ice thickness distribution initialization
319      !
320      CALL lim_hdf_init                ! set ice horizontal diffusion computation parameters
321      !
322      CALL lim_thd_init                ! set ice thermodynics parameters
323      !
324      CALL lim_thd_sal_init            ! set ice salinity parameters
325       
326      ! MV MP 2016
327      CALL lim_mp_init                 ! set melt ponds parameters
328      ! END MV MP 2016
329
330      IF( ln_limdyn )   CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation
331      !                                ! Initial sea-ice state
332      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
333         numit = 0
334         numit = nit000 - 1
335         CALL lim_istate
336      ELSE                                    ! start from a restart file
337         CALL lim_rst_read
338         numit = nit000 - 1
339      ENDIF
340      CALL lim_var_agg(2)
341      CALL lim_var_glo2eqv
342      !
343      CALL lim_sbc_init                 ! ice surface boundary condition
344      !
345      IF( ln_limdiahsb) CALL lim_diahsb_init  ! initialization for diags
346      !
347      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction
348      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
349      !
350      DO jj = 1, jpj
351         DO ji = 1, jpi
352            IF( gphit(ji,jj) > 0._wp ) THEN  ;  rn_amax_2d(ji,jj) = rn_amax_n  ! NH
353            ELSE                             ;  rn_amax_2d(ji,jj) = rn_amax_s  ! SH
354            ENDIF
355         END DO
356      END DO
357      !
358      nstart = numit  + nn_fsbc
359      nitrun = nitend - nit000 + 1
360      nlast  = numit  + nitrun
361      !
362      IF( nstock == 0 )   nstock = nlast + 1
363      !
364   END SUBROUTINE sbc_lim_init
365
366
367   SUBROUTINE lim_run_init
368      !!-------------------------------------------------------------------
369      !!                  ***  ROUTINE lim_run_init ***
370      !!
371      !! ** Purpose :   Definition some run parameter for ice model
372      !!
373      !! ** Method  :   Read the namicerun namelist and check the parameter
374      !!              values called at the first timestep (nit000)
375      !!
376      !! ** input   :   Namelist namicerun
377      !!-------------------------------------------------------------------
378      INTEGER  ::   ios                 ! Local integer output status for namelist read
379      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   &
380         &                cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice 
381      NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt 
382      !!-------------------------------------------------------------------
383      !
384      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice
385      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
386901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
387
388      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice
389      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
390902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
391      IF(lwm) WRITE ( numoni, namicerun )
392      !
393      REWIND( numnam_ice_ref )              ! Namelist namicediag in reference namelist : Parameters for ice
394      READ  ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
395903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
396
397      REWIND( numnam_ice_cfg )              ! Namelist namicediag in configuration namelist : Parameters for ice
398      READ  ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
399904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
400      IF(lwm) WRITE ( numoni, namicediag )
401      !
402      IF(lwp) THEN                        ! control print
403         WRITE(numout,*)
404         WRITE(numout,*) 'lim_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
405         WRITE(numout,*) ' ~~~~~~'
406         WRITE(numout,*) '   number of ice  categories                               = ', jpl
407         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i
408         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s
409         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n 
410         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s
411         WRITE(numout,*) '   Ice thermodynamics (T) or not (F)            ln_limthd  = ', ln_limthd
412         WRITE(numout,*) '   Ice dynamics       (T) or not (F)            ln_limdyn  = ', ln_limdyn
413         WRITE(numout,*) '     (ln_limdyn=T) Ice dynamics switch          nn_limdyn  = ', nn_limdyn
414         WRITE(numout,*) '       2: total'
415         WRITE(numout,*) '       1: advection only (no diffusion, no ridging/rafting)'
416         WRITE(numout,*) '       0: advection only (as 1 + prescribed velocity, bypass rheology)'
417         WRITE(numout,*) '     (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0)     = ', rn_uice
418         WRITE(numout,*) '     (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0)     = ', rn_vice
419         WRITE(numout,*)
420         WRITE(numout,*) '...and ice diagnostics'
421         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
422         WRITE(numout,*) '   Diagnose online heat/mass/salt budget     ln_limdiachk  = ', ln_limdiachk
423         WRITE(numout,*) '   Output          heat/mass/salt budget     ln_limdiahsb  = ', ln_limdiahsb
424         WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl
425         WRITE(numout,*) '   i-index for control prints (ln_limctl=true)             = ', iiceprt
426         WRITE(numout,*) '   j-index for control prints (ln_limctl=true)             = ', jiceprt
427      ENDIF
428      !
429      ! sea-ice timestep and inverse
430      rdt_ice   = REAL(nn_fsbc) * rdt 
431      r1_rdtice = 1._wp / rdt_ice 
432
433      ! inverse of nlay_i and nlay_s
434      r1_nlay_i = 1._wp / REAL( nlay_i, wp )
435      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
436      !
437      IF( lwp .AND. ln_bdy .AND. ln_limdiachk )  &
438         &   CALL ctl_warn('online conservation check activated but it does not work with BDY')
439      !
440      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice
441      !
442   END SUBROUTINE lim_run_init
443
444
445   SUBROUTINE lim_itd_init
446      !!------------------------------------------------------------------
447      !!                ***  ROUTINE lim_itd_init ***
448      !!
449      !! ** Purpose :   Initializes the ice thickness distribution
450      !! ** Method  :   ...
451      !! ** input   :   Namelist namiceitd
452      !!-------------------------------------------------------------------
453      INTEGER  ::   ios                 ! Local integer output status for namelist read
454      NAMELIST/namiceitd/ nn_catbnd, rn_himean
455      !
456      INTEGER  ::   jl                   ! dummy loop index
457      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
458      REAL(wp) ::   zhmax, znum, zden, zalpha !
459      !!------------------------------------------------------------------
460      !
461      REWIND( numnam_ice_ref )              ! Namelist namiceitd in reference namelist : Parameters for ice
462      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
463905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
464
465      REWIND( numnam_ice_cfg )              ! Namelist namiceitd in configuration namelist : Parameters for ice
466      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
467906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
468      IF(lwm) WRITE ( numoni, namiceitd )
469      !
470      IF(lwp) THEN                        ! control print
471         WRITE(numout,*)
472         WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution '
473         WRITE(numout,*) '~~~~~~~~~~~~'
474         WRITE(numout,*) '   shape of ice categories distribution                          nn_catbnd = ', nn_catbnd
475         WRITE(numout,*) '   mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean
476      ENDIF
477      !
478      !----------------------------------
479      !- Thickness categories boundaries
480      !----------------------------------
481      !
482      hi_max(:) = 0._wp
483      !
484      SELECT CASE ( nn_catbnd  )    ! type of ice categories distribution
485      !
486      CASE (1)                            !==  tanh function (CICE)  ==!
487         zc1 =  3._wp / REAL( jpl, wp )
488         zc2 = 10._wp * zc1
489         zc3 =  3._wp
490         DO jl = 1, jpl
491            zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
492            hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
493         END DO
494         !
495      CASE (2)                            !==  h^(-alpha) function  ==!
496         zalpha = 0.05_wp
497         zhmax  = 3._wp * rn_himean
498         DO jl = 1, jpl
499            znum = jpl * ( zhmax+1 )**zalpha
500            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
501            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
502         END DO
503         !
504      END SELECT
505      !
506      DO jl = 1, jpl                ! mean thickness by category
507         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
508      END DO
509      !
510      hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
511      !
512      IF(lwp) WRITE(numout,*)
513      IF(lwp) WRITE(numout,*) '      Thickness category boundaries '
514      IF(lwp) WRITE(numout,*) '         hi_max ', hi_max(0:jpl)
515      !
516   END SUBROUTINE lim_itd_init
517
518
519   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
520      !!---------------------------------------------------------------------
521      !!                  ***  ROUTINE ice_lim_flx  ***
522      !!
523      !! ** Purpose :   update the ice surface boundary condition by averaging and / or
524      !!                redistributing fluxes on ice categories
525      !!
526      !! ** Method  :   average then redistribute
527      !!
528      !! ** Action  :
529      !!---------------------------------------------------------------------
530      INTEGER                   , INTENT(in   ) ::   k_limflx   ! =-1 do nothing; =0 average ;
531      !                                                         ! = 1 average and redistribute ; =2 redistribute
532      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature
533      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
534      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux
535      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux
536      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity
537      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation
538      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
539      !
540      INTEGER  ::   jl      ! dummy loop index
541      !
542      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m    ! Mean albedo over all categories
543      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_m    ! Mean temperature over all categories
544      !
545      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m   ! Mean solar heat flux over all categories
546      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m   ! Mean non solar heat flux over all categories
547      REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m  ! Mean sublimation over all categories
548      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m   ! Mean d(qns)/dT over all categories
549      REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories
550      !!----------------------------------------------------------------------
551      !
552      IF( nn_timing == 1 )  CALL timing_start('ice_lim_flx')
553      !
554      SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==!
555      CASE( 0 , 1 )
556         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
557         !
558         z_qns_m  (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
559         z_qsr_m  (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
560         z_dqn_m  (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
561         z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
562         z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
563         DO jl = 1, jpl
564            pdqn_ice  (:,:,jl) = z_dqn_m(:,:)
565            pdevap_ice(:,:,jl) = z_devap_m(:,:)
566         END DO
567         !
568         DO jl = 1, jpl
569            pqns_ice (:,:,jl) = z_qns_m(:,:)
570            pqsr_ice (:,:,jl) = z_qsr_m(:,:)
571            pevap_ice(:,:,jl) = z_evap_m(:,:)
572         END DO
573         !
574         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
575      END SELECT
576      !
577      SELECT CASE( k_limflx )                              !==  redistribution on all ice categories  ==!
578      CASE( 1 , 2 )
579         CALL wrk_alloc( jpi,jpj, zalb_m, ztem_m )
580         !
581         zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )
582         ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) )
583         DO jl = 1, jpl
584            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
585            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
586            pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
587         END DO
588         !
589         CALL wrk_dealloc( jpi,jpj, zalb_m, ztem_m )
590      END SELECT
591      !
592      IF( nn_timing == 1 )  CALL timing_stop('ice_lim_flx')
593      !
594   END SUBROUTINE ice_lim_flx
595
596
597   SUBROUTINE sbc_lim_bef
598      !!----------------------------------------------------------------------
599      !!                  ***  ROUTINE sbc_lim_bef  ***
600      !!
601      !! ** purpose :  store ice variables at "before" time step
602      !!----------------------------------------------------------------------
603      a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area
604      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
605      v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume
606      v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume
607      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
608      smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content
609      oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content
610      u_ice_b(:,:)     = u_ice(:,:)
611      v_ice_b(:,:)     = v_ice(:,:)
612      !
613      at_i_b (:,:)     = SUM( a_i_b(:,:,:), dim=3 )
614     
615   END SUBROUTINE sbc_lim_bef
616
617
618   SUBROUTINE sbc_lim_diag0
619      !!----------------------------------------------------------------------
620      !!                  ***  ROUTINE sbc_lim_diag0  ***
621      !!
622      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
623      !!               of the time step
624      !!----------------------------------------------------------------------
625      sfx    (:,:) = 0._wp   ;
626      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp
627      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
628      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
629      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
630      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
631      !
632      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
633      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
634      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
635      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
636      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
637      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp 
638      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
639      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
640
641      ! MV MP 2016
642      wfx_pnd(:,:) = 0._wp   ;   wfx_snw_sum(:,:) = 0._wp
643      ! END MV MP 2016
644     
645      hfx_thd(:,:) = 0._wp   ;
646      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
647      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
648      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
649      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
650      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
651      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
652      hfx_err_dif(:,:) = 0._wp
653      wfx_err_sub(:,:) = 0._wp
654      !
655      afx_tot(:,:) = 0._wp   ;
656      afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp
657      !
658      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp
659      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
660
661      ! SIMIP diagnostics
662      diag_dms_dyn(:,:)  = 0._wp ; diag_dmi_dyn(:,:)  = 0._wp
663      diag_fc_bo(:,:)    = 0._wp ; diag_fc_su(:,:)    = 0._wp
664
665      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
666     
667   END SUBROUTINE sbc_lim_diag0
668
669
670   FUNCTION fice_cell_ave ( ptab )
671      !!--------------------------------------------------------------------------
672      !! * Compute average over categories, for grid cell (ice covered and free ocean)
673      !!--------------------------------------------------------------------------
674      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
675      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
676      INTEGER :: jl ! Dummy loop index
677
678      fice_cell_ave (:,:) = 0._wp
679      DO jl = 1, jpl
680         fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
681      END DO
682
683   END FUNCTION fice_cell_ave
684
685
686   FUNCTION fice_ice_ave ( ptab )
687      !!--------------------------------------------------------------------------
688      !! * Compute average over categories, for ice covered part of grid cell
689      !!--------------------------------------------------------------------------
690      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
691      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
692
693      fice_ice_ave (:,:) = 0.0_wp
694      WHERE ( at_i (:,:) > 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
695
696   END FUNCTION fice_ice_ave
697
698#else
699   !!----------------------------------------------------------------------
700   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
701   !!----------------------------------------------------------------------
702CONTAINS
703   !
704   SUBROUTINE sbc_ice_lim ( kt, ksbc )     ! Dummy routine
705      INTEGER, INTENT(in) ::   kt, ksbc
706      WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt
707   END SUBROUTINE sbc_ice_lim
708   !
709   SUBROUTINE sbc_lim_init                 ! Dummy routine
710   END SUBROUTINE sbc_lim_init
711#endif
712
713   !!======================================================================
714END MODULE sbcice_lim
Note: See TracBrowser for help on using the repository browser.