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.
sbcice_lim.F90 in branches/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 8809

Last change on this file since 8809 was 7777, checked in by clem, 7 years ago

fix for ticket #1765

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