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 @ 8360

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

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