New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limvar.F90 in branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

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

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