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/2015/dev_r5171_CNRS_LIM3_ice_atm/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2015/dev_r5171_CNRS_LIM3_ice_atm/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 5588

Last change on this file since 5588 was 5167, checked in by clem, 9 years ago

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