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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 8239

Last change on this file since 8239 was 8239, checked in by clem, 7 years ago

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