source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 6398

Last change on this file since 6398 was 6398, checked in by clem, 5 years ago

solve issue for out of bounds ice thicknesses. see ticket #1695

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