source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 8142

Last change on this file since 8142 was 8142, checked in by vancop, 3 years ago

Melt pond interfaces practically operational

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