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

Last change on this file since 8413 was 8413, checked in by clem, 3 years ago

continue changing names (again)

  • Property svn:keywords set to Id
File size: 34.4 KB
Line 
1MODULE limvar
2   !!======================================================================
3   !!                       ***  MODULE limvar ***
4   !!                 Different sets of ice model variables
5   !!                   how to switch from one to another
6   !!
7   !!                 There are three sets of variables
8   !!                 VGLO : global variables of the model
9   !!                        - v_i (jpi,jpj,jpl)
10   !!                        - v_s (jpi,jpj,jpl)
11   !!                        - a_i (jpi,jpj,jpl)
12   !!                        - t_s (jpi,jpj,jpl)
13   !!                        - e_i (jpi,jpj,nlay_i,jpl)
14   !!                        - smv_i(jpi,jpj,jpl)
15   !!                        - oa_i (jpi,jpj,jpl)
16   !!                 VEQV : equivalent variables sometimes used in the model
17   !!                        - ht_i(jpi,jpj,jpl)
18   !!                        - ht_s(jpi,jpj,jpl)
19   !!                        - t_i (jpi,jpj,nlay_i,jpl)
20   !!                        ...
21   !!                 VAGG : aggregate variables, averaged/summed over all
22   !!                        thickness categories
23   !!                        - vt_i(jpi,jpj)
24   !!                        - vt_s(jpi,jpj)
25   !!                        - at_i(jpi,jpj)
26   !!                        - et_s(jpi,jpj)  !total snow heat content
27   !!                        - et_i(jpi,jpj)  !total ice thermal content
28   !!                        - smt_i(jpi,jpj) !mean ice salinity
29   !!                        - tm_i (jpi,jpj) !mean ice temperature
30   !!======================================================================
31   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
32   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
33   !!----------------------------------------------------------------------
34#if defined key_lim3
35   !!----------------------------------------------------------------------
36   !!   'key_lim3'                                      LIM3 sea-ice model
37   !!----------------------------------------------------------------------
38   USE par_oce        ! ocean parameters
39   USE phycst         ! physical constants (ocean directory)
40   USE sbc_oce , ONLY : sss_m
41   USE ice            ! ice variables
42   USE thd_ice        ! ice variables (thermodynamics)
43   !
44   USE in_out_manager ! I/O manager
45   USE lib_mpp        ! MPP library
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), DIMENSION(jpi,jpj,jpl) ::   z_slope_s, zalpha
305      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
306      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
307      !!------------------------------------------------------------------
308
309      !---------------------------------------
310      ! Vertically constant, constant in time
311      !---------------------------------------
312      IF(  nn_icesal == 1  )  THEN
313         s_i (:,:,:,:) = rn_icesal
314         sm_i(:,:,:)   = rn_icesal
315      ENDIF
316
317      !-----------------------------------
318      ! Salinity profile, varying in time
319      !-----------------------------------
320      IF(  nn_icesal == 2  ) THEN
321         !
322         DO jk = 1, nlay_i
323            s_i(:,:,jk,:)  = sm_i(:,:,:)
324         END DO
325         !
326         DO jl = 1, jpl                               ! Slope of the linear profile
327            DO jj = 1, jpj
328               DO ji = 1, jpi
329                  rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) )
330                  z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) )
331               END DO
332            END DO
333         END DO
334         !
335         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
336         zfac1 = zsi1  / ( zsi1 - zsi0 )
337         !
338         zalpha(:,:,:) = 0._wp
339         DO jl = 1, jpl
340            DO jj = 1, jpj
341               DO ji = 1, jpi
342                  ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
343                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , zsi0 - sm_i(ji,jj,jl) ) ) 
344                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
345                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) ) 
346                  ! If 2.sm_i GE sss_m then rswitch = 1
347                  ! this is to force a constant salinity profile in the Baltic Sea
348                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
349                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
350                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
351               END DO
352            END DO
353         END DO
354
355         ! Computation of the profile
356         DO jl = 1, jpl
357            DO jk = 1, nlay_i
358               DO jj = 1, jpj
359                  DO ji = 1, jpi
360                     !                                      ! linear profile with 0 at the surface
361                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
362                     !                                      ! weighting the profile
363                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
364                     !                                      ! bounding salinity
365                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) )
366                  END DO
367               END DO
368            END DO
369         END DO
370         !
371      ENDIF ! nn_icesal
372
373      !-------------------------------------------------------
374      ! Vertically varying salinity profile, constant in time
375      !-------------------------------------------------------
376
377      IF(  nn_icesal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
378         !
379         sm_i(:,:,:) = 2.30_wp
380         !
381         DO jl = 1, jpl
382            DO jk = 1, nlay_i
383               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
384               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
385               s_i(:,:,jk,jl) =  zsal
386            END DO
387         END DO
388         !
389      ENDIF ! nn_icesal
390      !
391   END SUBROUTINE lim_var_salprof
392
393
394   SUBROUTINE lim_var_bv
395      !!------------------------------------------------------------------
396      !!                ***  ROUTINE lim_var_bv ***
397      !!
398      !! ** Purpose :   computes mean brine volume (%) in sea ice
399      !!
400      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
401      !!
402      !! References : Vancoppenolle et al., JGR, 2007
403      !!------------------------------------------------------------------
404      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
405      !!------------------------------------------------------------------
406      !
407      bvm_i(:,:)   = 0._wp
408      bv_i (:,:,:) = 0._wp
409      DO jl = 1, jpl
410         DO jk = 1, nlay_i
411            DO jj = 1, jpj
412               DO ji = 1, jpi
413                  rswitch        = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
414                  bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i  &
415                     &                            / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
416               END DO
417            END DO
418         END DO
419         
420         DO jj = 1, jpj
421            DO ji = 1, jpi
422               rswitch      = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
423               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 )
424            END DO
425         END DO
426      END DO
427      !
428   END SUBROUTINE lim_var_bv
429
430
431   SUBROUTINE lim_var_salprof1d
432      !!-------------------------------------------------------------------
433      !!                  ***  ROUTINE lim_thd_salprof1d  ***
434      !!
435      !! ** Purpose :   1d computation of the sea ice salinity profile
436      !!                Works with 1d vectors and is used by thermodynamic modules
437      !!-------------------------------------------------------------------
438      INTEGER  ::   ji, jk    ! dummy loop indices
439      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars
440      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      -
441      !
442      REAL(wp), DIMENSION(jpij) ::   z_slope_s
443      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
444      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
445      !!---------------------------------------------------------------------
446
447      !---------------------------------------
448      ! Vertically constant, constant in time
449      !---------------------------------------
450      IF( nn_icesal == 1 )   s_i_1d(:,:) = rn_icesal
451
452      !------------------------------------------------------
453      ! Vertically varying salinity profile, varying in time
454      !------------------------------------------------------
455
456      IF(  nn_icesal == 2  ) THEN
457         !
458         DO ji = 1, nidx          ! Slope of the linear profile zs_zero
459            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
460            z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
461         END DO
462
463         ! Weighting factor between zs_zero and zs_inf
464         !---------------------------------------------
465         zfac0 = 1._wp / ( zsi0 - zsi1 )
466         zfac1 = zsi1 / ( zsi1 - zsi0 )
467         DO jk = 1, nlay_i
468            DO ji = 1, nidx
469               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
470               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
471               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
472               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
473               ! if 2.sm_i GE sss_m then rswitch = 1
474               ! this is to force a constant salinity profile in the Baltic Sea
475               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_1d(ji) ) )
476               !
477               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch )
478               !
479               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
480               ! weighting the profile
481               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
482               ! bounding salinity
483               s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) )
484            END DO
485         END DO
486
487      ENDIF 
488
489      !-------------------------------------------------------
490      ! Vertically varying salinity profile, constant in time
491      !-------------------------------------------------------
492
493      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
494         !
495         sm_i_1d(:) = 2.30_wp
496         !
497         DO jk = 1, nlay_i
498            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
499            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
500            DO ji = 1, nidx
501               s_i_1d(ji,jk) = zsal
502            END DO
503         END DO
504         !
505      ENDIF
506      !
507   END SUBROUTINE lim_var_salprof1d
508
509   SUBROUTINE lim_var_zapsmall
510      !!-------------------------------------------------------------------
511      !!                   ***  ROUTINE lim_var_zapsmall ***
512      !!
513      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
514      !!
515      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
516      !!-------------------------------------------------------------------
517      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
518      REAL(wp) ::   zsal, zvi, zvs, zei, zes, zvp
519      !!-------------------------------------------------------------------
520      DO jl = 1, jpl
521
522         !-----------------------------------------------------------------
523         ! Zap ice energy and use ocean heat to melt ice
524         !-----------------------------------------------------------------
525         DO jk = 1, nlay_i
526            DO jj = 1 , jpj
527               DO ji = 1 , jpi
528                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
529                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
530                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
531                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
532                  zei              = e_i(ji,jj,jk,jl)
533                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
534                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
535                  ! update exchanges with ocean
536                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
537               END DO
538            END DO
539         END DO
540
541         DO jj = 1 , jpj
542            DO ji = 1 , jpi
543               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
544               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
545               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
546                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
547               zsal = smv_i(ji,jj,  jl)
548               zvi  = v_i  (ji,jj,  jl)
549               zvs  = v_s  (ji,jj,  jl)
550               zes  = e_s  (ji,jj,1,jl)
551               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw )   zvp  = v_ip (ji,jj  ,jl)
552               !-----------------------------------------------------------------
553               ! Zap snow energy
554               !-----------------------------------------------------------------
555               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
556               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
557
558               !-----------------------------------------------------------------
559               ! zap ice and snow volume, add water and salt to ocean
560               !-----------------------------------------------------------------
561               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
562               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
563               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
564               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
565               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
566               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
567               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
568
569               ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch
570               ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch
571
572               ! MV MP 2016
573               IF ( ln_pnd ) THEN
574                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
575                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
576                  IF ( ln_pnd_fw )   wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice
577               ENDIF
578               ! END MV MP 2016
579
580               ! update exchanges with ocean
581               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
582               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
583               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
584               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
585            END DO
586         END DO
587      END DO 
588
589      ! to be sure that at_i is the sum of a_i(jl)
590      at_i (:,:) = 0._wp
591      DO jl = 1, jpl
592         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
593      END DO
594
595      ! open water = 1 if at_i=0
596      DO jj = 1, jpj
597         DO ji = 1, jpi
598            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
599            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
600         END DO
601      END DO
602
603      !
604   END SUBROUTINE lim_var_zapsmall
605
606   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
607      !!------------------------------------------------------------------
608      !!                ***  ROUTINE lim_var_itd   ***
609      !!
610      !! ** Purpose :  converting 1-cat ice to multiple ice categories
611      !!
612      !!                  ice thickness distribution follows a gaussian law
613      !!               around the concentration of the most likely ice thickness
614      !!                           (similar as iceistate.F90)
615      !!
616      !! ** Method:   Iterative procedure
617      !!               
618      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
619      !!
620      !!               2) Check whether the distribution conserves area and volume, positivity and
621      !!                  category boundaries
622      !!             
623      !!               3) If not (input ice is too thin), the last category is empty and
624      !!                  the number of categories is reduced (jpl-1)
625      !!
626      !!               4) Iterate until ok (SUM(itest(:) = 4)
627      !!
628      !! ** Arguments : zhti: 1-cat ice thickness
629      !!                zhts: 1-cat snow depth
630      !!                zai : 1-cat ice concentration
631      !!
632      !! ** Output    : jpl-cat
633      !!
634      !!  (Example of application: BDY forcings when input are cell averaged) 
635      !!
636      !!-------------------------------------------------------------------
637      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
638      !!                    2014    (C. Rousset)        Rewriting
639      !!-------------------------------------------------------------------
640      !! Local variables
641      INTEGER  :: ji, jk, jl             ! dummy loop indices
642      INTEGER  :: ijpij, i_fill, jl0 
643      REAL(wp) :: zarg, zV, zconv, zdh, zdv
644      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
645      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
646      INTEGER , DIMENSION(4)                  ::   itest
647 
648      !--------------------------------------------------------------------
649      ! initialisation of variables
650      !--------------------------------------------------------------------
651      ijpij = SIZE(zhti,1)
652      zht_i(1:ijpij,1:jpl) = 0._wp
653      zht_s(1:ijpij,1:jpl) = 0._wp
654      za_i (1:ijpij,1:jpl) = 0._wp
655
656      ! ----------------------------------------
657      ! distribution over the jpl ice categories
658      ! ----------------------------------------
659      DO ji = 1, ijpij
660         
661         IF( zhti(ji) > 0._wp ) THEN
662
663            ! find which category (jl0) the input ice thickness falls into
664            jl0 = jpl
665            DO jl = 1, jpl
666               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
667                  jl0 = jl
668                  CYCLE
669               ENDIF
670            END DO
671
672            ! initialisation of tests
673            itest(:)  = 0
674         
675            i_fill = jpl + 1                                             !====================================
676            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
677               ! iteration                                               !====================================
678               i_fill = i_fill - 1
679               
680               ! initialisation of ice variables for each try
681               zht_i(ji,1:jpl) = 0._wp
682               za_i (ji,1:jpl) = 0._wp
683               itest(:)        = 0     
684               
685               ! *** case very thin ice: fill only category 1
686               IF ( i_fill == 1 ) THEN
687                  zht_i(ji,1) = zhti(ji)
688                  za_i (ji,1) = zai (ji)
689                 
690               ! *** case ice is thicker: fill categories >1
691               ELSE
692
693                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean
694                  DO jl = 1, i_fill - 1
695                     zht_i(ji,jl) = hi_mean(jl)
696                  END DO
697                 
698                  ! Concentrations in the (i_fill-1) categories
699                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
700                  DO jl = 1, i_fill - 1
701                     IF ( jl /= jl0 ) THEN
702                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
703                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
704                     ENDIF
705                  END DO
706                 
707                  ! Concentration in the last (i_fill) category
708                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
709                 
710                  ! Ice thickness in the last (i_fill) category
711                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
712                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
713                 
714                  ! clem: correction if concentration of upper cat is greater than lower cat
715                  !       (it should be a gaussian around jl0 but sometimes it is not)
716                  IF ( jl0 /= jpl ) THEN
717                     DO jl = jpl, jl0+1, -1
718                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
719                           zdv = zht_i(ji,jl) * za_i(ji,jl)
720                           zht_i(ji,jl    ) = 0._wp
721                           za_i (ji,jl    ) = 0._wp
722                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
723                        END IF
724                     ENDDO
725                  ENDIF
726               
727               ENDIF ! case ice is thick or thin
728           
729               !---------------------
730               ! Compatibility tests
731               !---------------------
732               ! Test 1: area conservation
733               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
734               IF ( zconv < epsi06 ) itest(1) = 1
735           
736               ! Test 2: volume conservation
737               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
738               IF ( zconv < epsi06 ) itest(2) = 1
739               
740               ! Test 3: thickness of the last category is in-bounds ?
741               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
742               
743               ! Test 4: positivity of ice concentrations
744               itest(4) = 1
745               DO jl = 1, i_fill
746                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
747               END DO
748               !                                         !============================
749            END DO                                       ! end iteration on categories
750               !                                         !============================
751         ENDIF ! if zhti > 0
752      END DO ! i loop
753
754      ! ------------------------------------------------
755      ! Adding Snow in each category where za_i is not 0
756      ! ------------------------------------------------
757      DO jl = 1, jpl
758         DO ji = 1, ijpij
759            IF( za_i(ji,jl) > 0._wp ) THEN
760               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
761               ! In case snow load is in excess that would lead to transformation from snow to ice
762               ! Then, transfer the snow excess into the ice (different from limthd_dh)
763               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
764               ! recompute ht_i, ht_s avoiding out of bounds values
765               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
766               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
767            ENDIF
768         ENDDO
769      ENDDO
770      !
771    END SUBROUTINE lim_var_itd
772
773
774#else
775   !!----------------------------------------------------------------------
776   !!   Default option         Dummy module          NO  LIM3 sea-ice model
777   !!----------------------------------------------------------------------
778CONTAINS
779   SUBROUTINE lim_var_agg          ! Empty routines
780   END SUBROUTINE lim_var_agg
781   SUBROUTINE lim_var_glo2eqv      ! Empty routines
782   END SUBROUTINE lim_var_glo2eqv
783   SUBROUTINE lim_var_eqv2glo      ! Empty routines
784   END SUBROUTINE lim_var_eqv2glo
785   SUBROUTINE lim_var_salprof      ! Empty routines
786   END SUBROUTINE lim_var_salprof
787   SUBROUTINE lim_var_bv           ! Emtpy routines
788   END SUBROUTINE lim_var_bv
789   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
790   END SUBROUTINE lim_var_salprof1d
791   SUBROUTINE lim_var_zapsmall
792   END SUBROUTINE lim_var_zapsmall
793   SUBROUTINE lim_var_itd
794   END SUBROUTINE lim_var_itd
795#endif
796
797   !!======================================================================
798END MODULE limvar
Note: See TracBrowser for help on using the repository browser.