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

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

solve issue #1712

  • Property svn:keywords set to Id
File size: 34.8 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  )  THEN
317         s_i (:,:,:,:) = rn_icesal
318         sm_i(:,:,:)   = rn_icesal
319      ENDIF
320
321      !-----------------------------------
322      ! Salinity profile, varying in time
323      !-----------------------------------
324      IF(  nn_icesal == 2  ) THEN
325         !
326         DO jk = 1, nlay_i
327            s_i(:,:,jk,:)  = sm_i(:,:,:)
328         END DO
329         !
330         DO jl = 1, jpl                               ! Slope of the linear profile
331            DO jj = 1, jpj
332               DO ji = 1, jpi
333                  rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) )
334                  z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) )
335               END DO
336            END DO
337         END DO
338         !
339         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
340         zfac1 = zsi1  / ( zsi1 - zsi0 )
341         !
342         zalpha(:,:,:) = 0._wp
343         DO jl = 1, jpl
344            DO jj = 1, jpj
345               DO ji = 1, jpi
346                  ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
347                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , zsi0 - sm_i(ji,jj,jl) ) ) 
348                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
349                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) ) 
350                  ! If 2.sm_i GE sss_m then rswitch = 1
351                  ! this is to force a constant salinity profile in the Baltic Sea
352                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
353                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
354                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
355               END DO
356            END DO
357         END DO
358
359         ! Computation of the profile
360         DO jl = 1, jpl
361            DO jk = 1, nlay_i
362               DO jj = 1, jpj
363                  DO ji = 1, jpi
364                     !                                      ! linear profile with 0 at the surface
365                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
366                     !                                      ! weighting the profile
367                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
368                     !                                      ! bounding salinity
369                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) )
370                  END DO
371               END DO
372            END DO
373         END DO
374         !
375      ENDIF ! nn_icesal
376
377      !-------------------------------------------------------
378      ! Vertically varying salinity profile, constant in time
379      !-------------------------------------------------------
380
381      IF(  nn_icesal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
382         !
383         sm_i(:,:,:) = 2.30_wp
384         !
385         DO jl = 1, jpl
386            DO jk = 1, nlay_i
387               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
388               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
389               s_i(:,:,jk,jl) =  zsal
390            END DO
391         END DO
392         !
393      ENDIF ! nn_icesal
394      !
395      CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
396      !
397   END SUBROUTINE lim_var_salprof
398
399
400   SUBROUTINE lim_var_icetm
401      !!------------------------------------------------------------------
402      !!                ***  ROUTINE lim_var_icetm ***
403      !!
404      !! ** Purpose :   computes mean sea ice temperature
405      !!------------------------------------------------------------------
406      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
407      !!------------------------------------------------------------------
408
409      ! Mean sea ice temperature
410      vt_i (:,:) = 0._wp
411      DO jl = 1, jpl
412         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
413      END DO
414
415      tm_i(:,:) = 0._wp
416      DO jl = 1, jpl
417         DO jk = 1, nlay_i
418            DO jj = 1, jpj
419               DO ji = 1, jpi
420                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
421                  tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
422                     &            / MAX( vt_i(ji,jj) , epsi10 )
423               END DO
424            END DO
425         END DO
426      END DO
427      tm_i = tm_i + rt0
428
429   END SUBROUTINE lim_var_icetm
430
431
432   SUBROUTINE lim_var_bv
433      !!------------------------------------------------------------------
434      !!                ***  ROUTINE lim_var_bv ***
435      !!
436      !! ** Purpose :   computes mean brine volume (%) in sea ice
437      !!
438      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
439      !!
440      !! References : Vancoppenolle et al., JGR, 2007
441      !!------------------------------------------------------------------
442      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
443      REAL(wp) ::   zbvi             ! local scalars
444      !!------------------------------------------------------------------
445      !
446      vt_i (:,:) = 0._wp
447      DO jl = 1, jpl
448         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
449      END DO
450
451      bv_i(:,:) = 0._wp
452      DO jl = 1, jpl
453         DO jk = 1, nlay_i
454            DO jj = 1, jpj
455               DO ji = 1, jpi
456                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
457                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   &
458                     &                   * v_i(ji,jj,jl) * r1_nlay_i
459                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) )  )
460                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi20 )
461               END DO
462            END DO
463         END DO
464      END DO
465      !
466   END SUBROUTINE lim_var_bv
467
468
469   SUBROUTINE lim_var_salprof1d( kideb, kiut )
470      !!-------------------------------------------------------------------
471      !!                  ***  ROUTINE lim_thd_salprof1d  ***
472      !!
473      !! ** Purpose :   1d computation of the sea ice salinity profile
474      !!                Works with 1d vectors and is used by thermodynamic modules
475      !!-------------------------------------------------------------------
476      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
477      !
478      INTEGER  ::   ji, jk    ! dummy loop indices
479      INTEGER  ::   ii, ij    ! local integers
480      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars
481      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      -
482      !
483      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
484      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
485      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
486      !!---------------------------------------------------------------------
487
488      CALL wrk_alloc( jpij, z_slope_s )
489
490      !---------------------------------------
491      ! Vertically constant, constant in time
492      !---------------------------------------
493      IF( nn_icesal == 1 )   s_i_1d(:,:) = rn_icesal
494
495      !------------------------------------------------------
496      ! Vertically varying salinity profile, varying in time
497      !------------------------------------------------------
498
499      IF(  nn_icesal == 2  ) THEN
500         !
501         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
502            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
503            z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
504         END DO
505
506         ! Weighting factor between zs_zero and zs_inf
507         !---------------------------------------------
508         zfac0 = 1._wp / ( zsi0 - zsi1 )
509         zfac1 = zsi1 / ( zsi1 - zsi0 )
510         DO jk = 1, nlay_i
511            DO ji = kideb, kiut
512               ii =  MOD( npb(ji) - 1 , jpi ) + 1
513               ij =     ( npb(ji) - 1 ) / jpi + 1
514               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
515               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
516               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
517               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
518               ! if 2.sm_i GE sss_m then rswitch = 1
519               ! this is to force a constant salinity profile in the Baltic Sea
520               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
521               !
522               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch )
523               !
524               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
525               ! weighting the profile
526               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
527               ! bounding salinity
528               s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) )
529            END DO
530         END DO
531
532      ENDIF 
533
534      !-------------------------------------------------------
535      ! Vertically varying salinity profile, constant in time
536      !-------------------------------------------------------
537
538      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
539         !
540         sm_i_1d(:) = 2.30_wp
541         !
542         DO jk = 1, nlay_i
543            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
544            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
545            DO ji = kideb, kiut
546               s_i_1d(ji,jk) = zsal
547            END DO
548         END DO
549         !
550      ENDIF
551      !
552      CALL wrk_dealloc( jpij, z_slope_s )
553      !
554   END SUBROUTINE lim_var_salprof1d
555
556   SUBROUTINE lim_var_zapsmall
557      !!-------------------------------------------------------------------
558      !!                   ***  ROUTINE lim_var_zapsmall ***
559      !!
560      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
561      !!
562      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
563      !!-------------------------------------------------------------------
564      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
565      REAL(wp) ::   zsal, zvi, zvs, zei, zes
566      !!-------------------------------------------------------------------
567      at_i (:,:) = 0._wp
568      DO jl = 1, jpl
569         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
570      END DO
571
572      DO jl = 1, jpl
573
574         !-----------------------------------------------------------------
575         ! Zap ice energy and use ocean heat to melt ice
576         !-----------------------------------------------------------------
577         DO jk = 1, nlay_i
578            DO jj = 1 , jpj
579               DO ji = 1 , jpi
580                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
581                  rswitch          = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
582                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
583                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
584                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
585                  zei              = e_i(ji,jj,jk,jl)
586                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
587                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
588                  ! update exchanges with ocean
589                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
590               END DO
591            END DO
592         END DO
593
594         DO jj = 1 , jpj
595            DO ji = 1 , jpi
596               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
597               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
598               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
599               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
600                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
601               zsal = smv_i(ji,jj,  jl)
602               zvi  = v_i  (ji,jj,  jl)
603               zvs  = v_s  (ji,jj,  jl)
604               zes  = e_s  (ji,jj,1,jl)
605               !-----------------------------------------------------------------
606               ! Zap snow energy
607               !-----------------------------------------------------------------
608               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
609               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
610
611               !-----------------------------------------------------------------
612               ! zap ice and snow volume, add water and salt to ocean
613               !-----------------------------------------------------------------
614               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
615               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
616               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
617               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
618               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
619               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
620               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
621
622               ! update exchanges with ocean
623               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
624               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
625               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
626               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
627            END DO
628         END DO
629      END DO 
630
631      ! to be sure that at_i is the sum of a_i(jl)
632      at_i (:,:) = 0._wp
633      DO jl = 1, jpl
634         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
635      END DO
636
637      ! open water = 1 if at_i=0
638      DO jj = 1, jpj
639         DO ji = 1, jpi
640            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
641            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
642         END DO
643      END DO
644
645      !
646   END SUBROUTINE lim_var_zapsmall
647
648   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
649      !!------------------------------------------------------------------
650      !!                ***  ROUTINE lim_var_itd   ***
651      !!
652      !! ** Purpose :  converting 1-cat ice to multiple ice categories
653      !!
654      !!                  ice thickness distribution follows a gaussian law
655      !!               around the concentration of the most likely ice thickness
656      !!                           (similar as limistate.F90)
657      !!
658      !! ** Method:   Iterative procedure
659      !!               
660      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
661      !!
662      !!               2) Check whether the distribution conserves area and volume, positivity and
663      !!                  category boundaries
664      !!             
665      !!               3) If not (input ice is too thin), the last category is empty and
666      !!                  the number of categories is reduced (jpl-1)
667      !!
668      !!               4) Iterate until ok (SUM(itest(:) = 4)
669      !!
670      !! ** Arguments : zhti: 1-cat ice thickness
671      !!                zhts: 1-cat snow depth
672      !!                zai : 1-cat ice concentration
673      !!
674      !! ** Output    : jpl-cat
675      !!
676      !!  (Example of application: BDY forcings when input are cell averaged) 
677      !!
678      !!-------------------------------------------------------------------
679      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
680      !!                    2014    (C. Rousset)        Rewriting
681      !!-------------------------------------------------------------------
682      !! Local variables
683      INTEGER  :: ji, jk, jl             ! dummy loop indices
684      INTEGER  :: ijpij, i_fill, jl0 
685      REAL(wp) :: zarg, zV, zconv, zdh
686      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
687      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
688      INTEGER , POINTER, DIMENSION(:)         ::   itest
689 
690      CALL wrk_alloc( 4, itest )
691      !--------------------------------------------------------------------
692      ! initialisation of variables
693      !--------------------------------------------------------------------
694      ijpij = SIZE(zhti,1)
695      zht_i(1:ijpij,1:jpl) = 0._wp
696      zht_s(1:ijpij,1:jpl) = 0._wp
697      za_i (1:ijpij,1:jpl) = 0._wp
698
699      ! ----------------------------------------
700      ! distribution over the jpl ice categories
701      ! ----------------------------------------
702      DO ji = 1, ijpij
703         
704         IF( zhti(ji) > 0._wp ) THEN
705
706         ! initialisation of tests
707         itest(:)  = 0
708         
709         i_fill = jpl + 1                                             !====================================
710         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
711            ! iteration                                               !====================================
712            i_fill = i_fill - 1
713           
714            ! initialisation of ice variables for each try
715            zht_i(ji,1:jpl) = 0._wp
716            za_i (ji,1:jpl) = 0._wp
717           
718            ! *** case very thin ice: fill only category 1
719            IF ( i_fill == 1 ) THEN
720               zht_i(ji,1) = zhti(ji)
721               za_i (ji,1) = zai (ji)
722
723            ! *** case ice is thicker: fill categories >1
724            ELSE
725
726               ! Fill ice thicknesses except the last one (i_fill) by hmean
727               DO jl = 1, i_fill - 1
728                  zht_i(ji,jl) = hi_mean(jl)
729               END DO
730               
731               ! find which category (jl0) the input ice thickness falls into
732               jl0 = i_fill
733               DO jl = 1, i_fill
734                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
735                     jl0 = jl
736           CYCLE
737                  ENDIF
738               END DO
739               
740               ! Concentrations in the (i_fill-1) categories
741               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
742               DO jl = 1, i_fill - 1
743                  IF ( jl == jl0 ) CYCLE
744                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
745                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
746               END DO
747               
748               ! Concentration in the last (i_fill) category
749               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
750               
751               ! Ice thickness in the last (i_fill) category
752               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
753               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
754               
755            ENDIF ! case ice is thick or thin
756           
757            !---------------------
758            ! Compatibility tests
759            !---------------------
760            ! Test 1: area conservation
761            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
762            IF ( zconv < epsi06 ) itest(1) = 1
763           
764            ! Test 2: volume conservation
765            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
766            IF ( zconv < epsi06 ) itest(2) = 1
767           
768            ! Test 3: thickness of the last category is in-bounds ?
769            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
770           
771            ! Test 4: positivity of ice concentrations
772            itest(4) = 1
773            DO jl = 1, i_fill
774               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
775            END DO           
776                                                           !============================
777         END DO                                            ! end iteration on categories
778                                                           !============================
779         ENDIF ! if zhti > 0
780      END DO ! i loop
781
782      ! ------------------------------------------------
783      ! Adding Snow in each category where za_i is not 0
784      ! ------------------------------------------------
785      DO jl = 1, jpl
786         DO ji = 1, ijpij
787            IF( za_i(ji,jl) > 0._wp ) THEN
788               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
789               ! In case snow load is in excess that would lead to transformation from snow to ice
790               ! Then, transfer the snow excess into the ice (different from limthd_dh)
791               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
792               ! recompute ht_i, ht_s avoiding out of bounds values
793               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
794               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
795            ENDIF
796         ENDDO
797      ENDDO
798
799      CALL wrk_dealloc( 4, itest )
800      !
801    END SUBROUTINE lim_var_itd
802
803
804#else
805   !!----------------------------------------------------------------------
806   !!   Default option         Dummy module          NO  LIM3 sea-ice model
807   !!----------------------------------------------------------------------
808CONTAINS
809   SUBROUTINE lim_var_agg          ! Empty routines
810   END SUBROUTINE lim_var_agg
811   SUBROUTINE lim_var_glo2eqv      ! Empty routines
812   END SUBROUTINE lim_var_glo2eqv
813   SUBROUTINE lim_var_eqv2glo      ! Empty routines
814   END SUBROUTINE lim_var_eqv2glo
815   SUBROUTINE lim_var_salprof      ! Empty routines
816   END SUBROUTINE lim_var_salprof
817   SUBROUTINE lim_var_bv           ! Emtpy routines
818   END SUBROUTINE lim_var_bv
819   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
820   END SUBROUTINE lim_var_salprof1d
821   SUBROUTINE lim_var_zapsmall
822   END SUBROUTINE lim_var_zapsmall
823   SUBROUTINE lim_var_itd
824   END SUBROUTINE lim_var_itd
825#endif
826
827   !!======================================================================
828END MODULE limvar
Note: See TracBrowser for help on using the repository browser.