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.
limvar.F90 in branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 6853

Last change on this file since 6853 was 6853, checked in by clem, 8 years ago

correct bugs on LIM3 rheology and outputs. Make sure the ice thickness distribution (used in initialization and bdy boundary conditions) is a gaussian.
  • Property svn:keywords set to Id
File size: 34.3 KB
Line 
1MODULE limvar
2   !!======================================================================
3   !!                       ***  MODULE limvar ***
4   !!                 Different sets of ice model variables
5   !!                   how to switch from one to another
6   !!
7   !!                 There are three sets of variables
8   !!                 VGLO : global variables of the model
9   !!                        - v_i (jpi,jpj,jpl)
10   !!                        - v_s (jpi,jpj,jpl)
11   !!                        - a_i (jpi,jpj,jpl)
12   !!                        - t_s (jpi,jpj,jpl)
13   !!                        - e_i (jpi,jpj,nlay_i,jpl)
14   !!                        - smv_i(jpi,jpj,jpl)
15   !!                        - oa_i (jpi,jpj,jpl)
16   !!                 VEQV : equivalent variables sometimes used in the model
17   !!                        - ht_i(jpi,jpj,jpl)
18   !!                        - ht_s(jpi,jpj,jpl)
19   !!                        - t_i (jpi,jpj,nlay_i,jpl)
20   !!                        ...
21   !!                 VAGG : aggregate variables, averaged/summed over all
22   !!                        thickness categories
23   !!                        - vt_i(jpi,jpj)
24   !!                        - vt_s(jpi,jpj)
25   !!                        - at_i(jpi,jpj)
26   !!                        - et_s(jpi,jpj)  !total snow heat content
27   !!                        - et_i(jpi,jpj)  !total ice thermal content
28   !!                        - smt_i(jpi,jpj) !mean ice salinity
29   !!                        - tm_i (jpi,jpj) !mean ice temperature
30   !!======================================================================
31   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
32   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
33   !!----------------------------------------------------------------------
34#if defined key_lim3
35   !!----------------------------------------------------------------------
36   !!   'key_lim3'                                      LIM3 sea-ice model
37   !!----------------------------------------------------------------------
38   USE par_oce        ! ocean parameters
39   USE phycst         ! physical constants (ocean directory)
40   USE sbc_oce        ! Surface boundary condition: ocean fields
41   USE ice            ! ice variables
42   USE thd_ice        ! ice variables (thermodynamics)
43   USE in_out_manager ! I/O manager
44   USE lib_mpp        ! MPP library
45   USE wrk_nemo       ! work arrays
46   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
47
48   IMPLICIT NONE
49   PRIVATE
50
51   PUBLIC   lim_var_agg         
52   PUBLIC   lim_var_glo2eqv     
53   PUBLIC   lim_var_eqv2glo     
54   PUBLIC   lim_var_salprof     
55   PUBLIC   lim_var_bv           
56   PUBLIC   lim_var_salprof1d   
57   PUBLIC   lim_var_zapsmall
58   PUBLIC   lim_var_itd
59
60   !!----------------------------------------------------------------------
61   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE lim_var_agg( kn )
68      !!------------------------------------------------------------------
69      !!                ***  ROUTINE lim_var_agg  ***
70      !!
71      !! ** Purpose :   aggregates ice-thickness-category variables to all-ice variables
72      !!              i.e. it turns VGLO into VAGG
73      !! ** Method  :
74      !!
75      !! ** Arguments : n = 1, at_i vt_i only
76      !!                n = 2 everything
77      !!
78      !! note : you could add an argument when you need only at_i, vt_i
79      !!        and when you need everything
80      !!------------------------------------------------------------------
81      INTEGER, INTENT( in ) ::   kn     ! =1 at_i & vt only ; = what is needed
82      !
83      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
84      !!------------------------------------------------------------------
85
86      ! integrated values
87      vt_i (:,:) = SUM( v_i, dim=3 )
88      vt_s (:,:) = SUM( v_s, dim=3 )
89      at_i (:,:) = SUM( a_i, dim=3 )
90      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
91      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
92
93      ! open water fraction
94      DO jj = 1, jpj
95         DO ji = 1, jpi
96            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   
97         END DO
98      END DO
99
100      IF( kn > 1 ) THEN
101
102         ! mean ice/snow thickness
103         DO jj = 1, jpj
104            DO ji = 1, jpi
105               rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
106               htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
107               htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
108            ENDDO
109         ENDDO
110
111         ! mean temperature (K), salinity and age
112         smt_i(:,:) = 0._wp
113         tm_i(:,:)  = 0._wp
114         tm_su(:,:) = 0._wp
115         om_i (:,:) = 0._wp
116         DO jl = 1, jpl
117           
118            DO jj = 1, jpj
119               DO ji = 1, jpi
120                  rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
121                  tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
122                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi10 )
123               END DO
124            END DO
125           
126            DO jk = 1, nlay_i
127               DO jj = 1, jpj
128                  DO ji = 1, jpi
129                     rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
130                     tm_i(ji,jj)  = tm_i(ji,jj)  + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
131                        &            / MAX( vt_i(ji,jj) , epsi10 )
132                     smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl)  &
133                        &            / MAX( vt_i(ji,jj) , epsi10 )
134                  END DO
135               END DO
136            END DO
137         END DO
138         tm_i  = tm_i  + rt0
139         tm_su = tm_su + rt0
140         !
141      ENDIF
142      !
143   END SUBROUTINE lim_var_agg
144
145
146   SUBROUTINE lim_var_glo2eqv
147      !!------------------------------------------------------------------
148      !!                ***  ROUTINE lim_var_glo2eqv ***
149      !!
150      !! ** Purpose :   computes equivalent variables as function of global variables
151      !!              i.e. it turns VGLO into VEQV
152      !!------------------------------------------------------------------
153      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
154      REAL(wp) ::   zq_i, zaaa, zbbb, zccc, zdiscrim     ! local scalars
155      REAL(wp) ::   ztmelts, zq_s, zfac1, zfac2   !   -      -
156      !!------------------------------------------------------------------
157
158      !-------------------------------------------------------
159      ! Ice thickness, snow thickness, ice salinity, ice age
160      !-------------------------------------------------------
161      DO jl = 1, jpl
162         DO jj = 1, jpj
163            DO ji = 1, jpi
164               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
165               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
166            END DO
167         END DO
168      END DO
169      ! Force the upper limit of ht_i to always be < hi_max (99 m).
170      DO jj = 1, jpj
171         DO ji = 1, jpi
172            rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
173            ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
174            a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
175         END DO
176      END DO
177
178      DO jl = 1, jpl
179         DO jj = 1, jpj
180            DO ji = 1, jpi
181               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
182               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
183               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
184            END DO
185         END DO
186      END DO
187     
188      IF(  nn_icesal == 2  )THEN
189         DO jl = 1, jpl
190            DO jj = 1, jpj
191               DO ji = 1, jpi
192                  rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
193                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch
194                  !                                      ! bounding salinity
195                  sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin )
196               END DO
197            END DO
198         END DO
199      ENDIF
200
201      CALL lim_var_salprof      ! salinity profile
202
203      !-------------------
204      ! Ice temperatures
205      !-------------------
206      DO jl = 1, jpl
207         DO jk = 1, nlay_i
208            DO jj = 1, jpj
209               DO ji = 1, jpi
210                  !                                                              ! Energy of melting q(S,T) [J.m-3]
211                  rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
212                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 
213                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0              ! Ice layer melt temperature
214                  !
215                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
216                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus
217                  zccc       =  lfus * (ztmelts-rt0)
218                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
219                  t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )
220                  t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) )  ! -100 < t_i < ztmelts
221               END DO
222            END DO
223         END DO
224      END DO
225
226      !--------------------
227      ! Snow temperatures
228      !--------------------
229      zfac1 = 1._wp / ( rhosn * cpic )
230      zfac2 = lfus / cpic 
231      DO jl = 1, jpl
232         DO jk = 1, nlay_s
233            DO jj = 1, jpj
234               DO ji = 1, jpi
235                  !Energy of melting q(S,T) [J.m-3]
236                  rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
237                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)
238                  !
239                  t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 )
240                  t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) )     ! -100 < t_s < rt0
241               END DO
242            END DO
243         END DO
244      END DO
245
246      ! integrated values
247      vt_i (:,:) = SUM( v_i, dim=3 )
248      vt_s (:,:) = SUM( v_s, dim=3 )
249      at_i (:,:) = SUM( a_i, dim=3 )
250
251      !
252   END SUBROUTINE lim_var_glo2eqv
253
254
255   SUBROUTINE lim_var_eqv2glo
256      !!------------------------------------------------------------------
257      !!                ***  ROUTINE lim_var_eqv2glo ***
258      !!
259      !! ** Purpose :   computes global variables as function of equivalent variables
260      !!                i.e. it turns VEQV into VGLO
261      !! ** Method  :
262      !!
263      !! ** History :  (01-2006) Martin Vancoppenolle, UCL-ASTR
264      !!------------------------------------------------------------------
265      !
266      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
267      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
268      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
269      !
270   END SUBROUTINE lim_var_eqv2glo
271
272
273   SUBROUTINE lim_var_salprof
274      !!------------------------------------------------------------------
275      !!                ***  ROUTINE lim_var_salprof ***
276      !!
277      !! ** Purpose :   computes salinity profile in function of bulk salinity     
278      !!
279      !! ** Method  : If bulk salinity greater than zsi1,
280      !!              the profile is assumed to be constant (S_inf)
281      !!              If bulk salinity lower than zsi0,
282      !!              the profile is linear with 0 at the surface (S_zero)
283      !!              If it is between zsi0 and zsi1, it is a
284      !!              alpha-weighted linear combination of s_inf and s_zero
285      !!
286      !! ** References : Vancoppenolle et al., 2007
287      !!------------------------------------------------------------------
288      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
289      REAL(wp) ::   zfac0, zfac1, zsal
290      REAL(wp) ::   zswi0, zswi01, zargtemp , zs_zero   
291      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha
292      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
293      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
294      !!------------------------------------------------------------------
295
296      CALL wrk_alloc( jpi, jpj, jpl, z_slope_s, zalpha )
297
298      !---------------------------------------
299      ! Vertically constant, constant in time
300      !---------------------------------------
301      IF(  nn_icesal == 1  )  THEN
302         s_i (:,:,:,:) = rn_icesal
303         sm_i(:,:,:)   = rn_icesal
304      ENDIF
305         
306      !-----------------------------------
307      ! Salinity profile, varying in time
308      !-----------------------------------
309      IF(  nn_icesal == 2  ) THEN
310         !
311         DO jk = 1, nlay_i
312            s_i(:,:,jk,:)  = sm_i(:,:,:)
313         END DO
314         !
315         DO jl = 1, jpl                               ! Slope of the linear profile
316            DO jj = 1, jpj
317               DO ji = 1, jpi
318                  rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) )
319                  z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) )
320               END DO
321            END DO
322         END DO
323         !
324         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
325         zfac1 = zsi1  / ( zsi1 - zsi0 )
326         !
327         zalpha(:,:,:) = 0._wp
328         DO jl = 1, jpl
329            DO jj = 1, jpj
330               DO ji = 1, jpi
331                  ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
332                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , zsi0 - sm_i(ji,jj,jl) ) ) 
333                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
334                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) ) 
335                  ! If 2.sm_i GE sss_m then rswitch = 1
336                  ! this is to force a constant salinity profile in the Baltic Sea
337                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
338                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
339                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
340               END DO
341            END DO
342         END DO
343
344         ! Computation of the profile
345         DO jl = 1, jpl
346            DO jk = 1, nlay_i
347               DO jj = 1, jpj
348                  DO ji = 1, jpi
349                     !                                      ! linear profile with 0 at the surface
350                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
351                     !                                      ! weighting the profile
352                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
353                     !                                      ! bounding salinity
354                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) )
355                  END DO
356               END DO
357            END DO
358         END DO
359         !
360      ENDIF ! nn_icesal
361
362      !-------------------------------------------------------
363      ! Vertically varying salinity profile, constant in time
364      !-------------------------------------------------------
365
366      IF(  nn_icesal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
367         !
368         sm_i(:,:,:) = 2.30_wp
369         !
370         DO jl = 1, jpl
371            DO jk = 1, nlay_i
372               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
373               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
374               s_i(:,:,jk,jl) =  zsal
375            END DO
376         END DO
377         !
378      ENDIF ! nn_icesal
379      !
380      CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
381      !
382   END SUBROUTINE lim_var_salprof
383
384
385   SUBROUTINE lim_var_bv
386      !!------------------------------------------------------------------
387      !!                ***  ROUTINE lim_var_bv ***
388      !!
389      !! ** Purpose :   computes mean brine volume (%) in sea ice
390      !!
391      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
392      !!
393      !! References : Vancoppenolle et al., JGR, 2007
394      !!------------------------------------------------------------------
395      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
396      !!------------------------------------------------------------------
397      !
398      bvm_i(:,:)   = 0._wp
399      bv_i (:,:,:) = 0._wp
400      DO jl = 1, jpl
401         DO jk = 1, nlay_i
402            DO jj = 1, jpj
403               DO ji = 1, jpi
404                  rswitch        = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
405                  bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i  &
406                     &                            / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
407               END DO
408            END DO
409         END DO
410         
411         DO jj = 1, jpj
412            DO ji = 1, jpi
413               rswitch      = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
414               bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 )
415            END DO
416         END DO
417      END DO
418      !
419   END SUBROUTINE lim_var_bv
420
421
422   SUBROUTINE lim_var_salprof1d( kideb, kiut )
423      !!-------------------------------------------------------------------
424      !!                  ***  ROUTINE lim_thd_salprof1d  ***
425      !!
426      !! ** Purpose :   1d computation of the sea ice salinity profile
427      !!                Works with 1d vectors and is used by thermodynamic modules
428      !!-------------------------------------------------------------------
429      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
430      !
431      INTEGER  ::   ji, jk    ! dummy loop indices
432      INTEGER  ::   ii, ij    ! local integers
433      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars
434      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      -
435      !
436      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
437      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
438      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
439      !!---------------------------------------------------------------------
440
441      CALL wrk_alloc( jpij, z_slope_s )
442
443      !---------------------------------------
444      ! Vertically constant, constant in time
445      !---------------------------------------
446      IF( nn_icesal == 1 )   s_i_1d(:,:) = rn_icesal
447
448      !------------------------------------------------------
449      ! Vertically varying salinity profile, varying in time
450      !------------------------------------------------------
451
452      IF(  nn_icesal == 2  ) THEN
453         !
454         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
455            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
456            z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
457         END DO
458
459         ! Weighting factor between zs_zero and zs_inf
460         !---------------------------------------------
461         zfac0 = 1._wp / ( zsi0 - zsi1 )
462         zfac1 = zsi1 / ( zsi1 - zsi0 )
463         DO jk = 1, nlay_i
464            DO ji = kideb, kiut
465               ii =  MOD( npb(ji) - 1 , jpi ) + 1
466               ij =     ( npb(ji) - 1 ) / jpi + 1
467               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
468               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
469               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
470               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
471               ! if 2.sm_i GE sss_m then rswitch = 1
472               ! this is to force a constant salinity profile in the Baltic Sea
473               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
474               !
475               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch )
476               !
477               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
478               ! weighting the profile
479               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
480               ! bounding salinity
481               s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) )
482            END DO
483         END DO
484
485      ENDIF 
486
487      !-------------------------------------------------------
488      ! Vertically varying salinity profile, constant in time
489      !-------------------------------------------------------
490
491      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
492         !
493         sm_i_1d(:) = 2.30_wp
494         !
495         DO jk = 1, nlay_i
496            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
497            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
498            DO ji = kideb, kiut
499               s_i_1d(ji,jk) = zsal
500            END DO
501         END DO
502         !
503      ENDIF
504      !
505      CALL wrk_dealloc( jpij, z_slope_s )
506      !
507   END SUBROUTINE lim_var_salprof1d
508
509   SUBROUTINE lim_var_zapsmall
510      !!-------------------------------------------------------------------
511      !!                   ***  ROUTINE lim_var_zapsmall ***
512      !!
513      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
514      !!
515      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
516      !!-------------------------------------------------------------------
517      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
518      REAL(wp) ::   zsal, zvi, zvs, zei, zes
519      !!-------------------------------------------------------------------
520      at_i (:,:) = 0._wp
521      DO jl = 1, jpl
522         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
523      END DO
524
525      DO jl = 1, jpl
526
527         !-----------------------------------------------------------------
528         ! Zap ice energy and use ocean heat to melt ice
529         !-----------------------------------------------------------------
530         DO jk = 1, nlay_i
531            DO jj = 1 , jpj
532               DO ji = 1 , jpi
533                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
534                  rswitch          = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
535                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
536                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
537                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
538                  zei              = e_i(ji,jj,jk,jl)
539                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
540                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
541                  ! update exchanges with ocean
542                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
543               END DO
544            END DO
545         END DO
546
547         DO jj = 1 , jpj
548            DO ji = 1 , jpi
549               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
550               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
551               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
552               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
553                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
554               zsal = smv_i(ji,jj,  jl)
555               zvi  = v_i  (ji,jj,  jl)
556               zvs  = v_s  (ji,jj,  jl)
557               zes  = e_s  (ji,jj,1,jl)
558               !-----------------------------------------------------------------
559               ! Zap snow energy
560               !-----------------------------------------------------------------
561               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
562               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
563
564               !-----------------------------------------------------------------
565               ! zap ice and snow volume, add water and salt to ocean
566               !-----------------------------------------------------------------
567               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
568               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
569               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
570               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
571               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
572               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
573               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
574
575               ! update exchanges with ocean
576               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
577               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
578               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
579               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
580            END DO
581         END DO
582      END DO 
583
584      ! to be sure that at_i is the sum of a_i(jl)
585      at_i (:,:) = 0._wp
586      DO jl = 1, jpl
587         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
588      END DO
589
590      ! open water = 1 if at_i=0
591      DO jj = 1, jpj
592         DO ji = 1, jpi
593            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
594            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
595         END DO
596      END DO
597
598      !
599   END SUBROUTINE lim_var_zapsmall
600
601   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
602      !!------------------------------------------------------------------
603      !!                ***  ROUTINE lim_var_itd   ***
604      !!
605      !! ** Purpose :  converting 1-cat ice to multiple ice categories
606      !!
607      !!                  ice thickness distribution follows a gaussian law
608      !!               around the concentration of the most likely ice thickness
609      !!                           (similar as limistate.F90)
610      !!
611      !! ** Method:   Iterative procedure
612      !!               
613      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
614      !!
615      !!               2) Check whether the distribution conserves area and volume, positivity and
616      !!                  category boundaries
617      !!             
618      !!               3) If not (input ice is too thin), the last category is empty and
619      !!                  the number of categories is reduced (jpl-1)
620      !!
621      !!               4) Iterate until ok (SUM(itest(:) = 4)
622      !!
623      !! ** Arguments : zhti: 1-cat ice thickness
624      !!                zhts: 1-cat snow depth
625      !!                zai : 1-cat ice concentration
626      !!
627      !! ** Output    : jpl-cat
628      !!
629      !!  (Example of application: BDY forcings when input are cell averaged) 
630      !!
631      !!-------------------------------------------------------------------
632      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
633      !!                    2014    (C. Rousset)        Rewriting
634      !!-------------------------------------------------------------------
635      !! Local variables
636      INTEGER  :: ji, jk, jl             ! dummy loop indices
637      INTEGER  :: ijpij, i_fill, jl0 
638      REAL(wp) :: zarg, zV, zconv, zdh, zdv
639      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
640      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
641      INTEGER , POINTER, DIMENSION(:)         ::   itest
642 
643      Call wrk_alloc( 4, itest )
644      !--------------------------------------------------------------------
645      ! initialisation of variables
646      !--------------------------------------------------------------------
647      ijpij = SIZE(zhti,1)
648      zht_i(1:ijpij,1:jpl) = 0._wp
649      zht_s(1:ijpij,1:jpl) = 0._wp
650      za_i (1:ijpij,1:jpl) = 0._wp
651
652      ! ----------------------------------------
653      ! distribution over the jpl ice categories
654      ! ----------------------------------------
655      DO ji = 1, ijpij
656         
657         IF( zhti(ji) > 0._wp ) THEN
658
659            ! find which category (jl0) the input ice thickness falls into
660            jl0 = jpl
661            DO jl = 1, jpl
662               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
663                  jl0 = jl
664                  CYCLE
665               ENDIF
666            END DO
667
668            ! initialisation of tests
669            itest(:)  = 0
670         
671            i_fill = jpl + 1                                             !====================================
672            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
673               ! iteration                                               !====================================
674               i_fill = i_fill - 1
675               
676               ! initialisation of ice variables for each try
677               zht_i(ji,1:jpl) = 0._wp
678               za_i (ji,1:jpl) = 0._wp
679               itest(:)        = 0           
680               
681               ! *** case very thin ice: fill only category 1
682               IF ( i_fill == 1 ) THEN
683                  zht_i(ji,1) = zhti(ji)
684                  za_i (ji,1) = zai (ji)
685                 
686               ! *** case ice is thicker: fill categories >1
687               ELSE
688
689                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean
690                  DO jl = 1, i_fill - 1
691                     zht_i(ji,jl) = hi_mean(jl)
692                  END DO
693                 
694                  ! Concentrations in the (i_fill-1) categories
695                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
696                  DO jl = 1, i_fill - 1
697                     IF ( jl /= jl0 ) THEN
698                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
699                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
700                     ENDIF
701                  END DO
702                 
703                  ! Concentration in the last (i_fill) category
704                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
705                 
706                  ! Ice thickness in the last (i_fill) category
707                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
708                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
709                 
710                  ! clem: correction if concentration of upper cat is greater than lower cat
711                  !       (it should be a gaussian around jl0 but sometimes it is not)
712                  IF ( jl0 /= jpl ) THEN
713                     DO jl = jpl, jl0+1, -1
714                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
715                           zdv = zht_i(ji,jl) * za_i(ji,jl)
716                           zht_i(ji,jl    ) = 0._wp
717                           za_i (ji,jl    ) = 0._wp
718                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
719                        END IF
720                     ENDDO
721                  ENDIF
722               
723               ENDIF ! case ice is thick or thin
724               
725               !---------------------
726               ! Compatibility tests
727               !---------------------
728               ! Test 1: area conservation
729               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
730               IF ( zconv < epsi06 ) itest(1) = 1
731           
732               ! Test 2: volume conservation
733               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
734               IF ( zconv < epsi06 ) itest(2) = 1
735               
736               ! Test 3: thickness of the last category is in-bounds ?
737               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
738               
739               ! Test 4: positivity of ice concentrations
740               itest(4) = 1
741               DO jl = 1, i_fill
742                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
743               END DO
744               !                                         !============================
745            END DO                                       ! end iteration on categories
746               !                                         !============================
747         ENDIF ! if zhti > 0
748      END DO ! i loop
749     
750      ! ------------------------------------------------
751      ! Adding Snow in each category where za_i is not 0
752      ! ------------------------------------------------
753      DO jl = 1, jpl
754         DO ji = 1, ijpij
755            IF( za_i(ji,jl) > 0._wp ) THEN
756               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
757               ! In case snow load is in excess that would lead to transformation from snow to ice
758               ! Then, transfer the snow excess into the ice (different from limthd_dh)
759               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
760               ! recompute ht_i, ht_s avoiding out of bounds values
761               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
762               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
763            ENDIF
764         ENDDO
765      ENDDO
766
767      CALL wrk_dealloc( 4, itest )
768      !
769    END SUBROUTINE lim_var_itd
770
771
772#else
773   !!----------------------------------------------------------------------
774   !!   Default option         Dummy module          NO  LIM3 sea-ice model
775   !!----------------------------------------------------------------------
776CONTAINS
777   SUBROUTINE lim_var_agg          ! Empty routines
778   END SUBROUTINE lim_var_agg
779   SUBROUTINE lim_var_glo2eqv      ! Empty routines
780   END SUBROUTINE lim_var_glo2eqv
781   SUBROUTINE lim_var_eqv2glo      ! Empty routines
782   END SUBROUTINE lim_var_eqv2glo
783   SUBROUTINE lim_var_salprof      ! Empty routines
784   END SUBROUTINE lim_var_salprof
785   SUBROUTINE lim_var_bv           ! Emtpy routines
786   END SUBROUTINE lim_var_bv
787   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
788   END SUBROUTINE lim_var_salprof1d
789   SUBROUTINE lim_var_zapsmall
790   END SUBROUTINE lim_var_zapsmall
791   SUBROUTINE lim_var_itd
792   END SUBROUTINE lim_var_itd
793#endif
794
795   !!======================================================================
796END MODULE limvar
Note: See TracBrowser for help on using the repository browser.