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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 5134

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

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