source: NEMO/releases/release-3.6/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 10379

Last change on this file since 10379 was 7814, checked in by clem, 4 years ago

cosmetic changes for better phasing with trunk

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