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

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

STEP3 (1): clean separation between sea-ice and ocean

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