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

source: branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 7067

Last change on this file since 7067 was 7067, checked in by cetlod, 8 years ago

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