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/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

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

Last change on this file since 8085 was 8085, checked in by vancop, 8 years ago

More on melt ponds in LIM

  • Property svn:keywords set to Id
File size: 34.9 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 ( nn_pnd_scheme > 0 ) 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               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               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
588               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
589               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice
590               ! END MV MP 2016
591
592               ! update exchanges with ocean
593               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
594               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
595               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
596               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
597            END DO
598         END DO
599      END DO 
600
601      ! to be sure that at_i is the sum of a_i(jl)
602      at_i (:,:) = 0._wp
603      DO jl = 1, jpl
604         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
605      END DO
606
607      ! open water = 1 if at_i=0
608      DO jj = 1, jpj
609         DO ji = 1, jpi
610            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
611            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
612         END DO
613      END DO
614
615      !
616   END SUBROUTINE lim_var_zapsmall
617
618   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
619      !!------------------------------------------------------------------
620      !!                ***  ROUTINE lim_var_itd   ***
621      !!
622      !! ** Purpose :  converting 1-cat ice to multiple ice categories
623      !!
624      !!                  ice thickness distribution follows a gaussian law
625      !!               around the concentration of the most likely ice thickness
626      !!                           (similar as limistate.F90)
627      !!
628      !! ** Method:   Iterative procedure
629      !!               
630      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
631      !!
632      !!               2) Check whether the distribution conserves area and volume, positivity and
633      !!                  category boundaries
634      !!             
635      !!               3) If not (input ice is too thin), the last category is empty and
636      !!                  the number of categories is reduced (jpl-1)
637      !!
638      !!               4) Iterate until ok (SUM(itest(:) = 4)
639      !!
640      !! ** Arguments : zhti: 1-cat ice thickness
641      !!                zhts: 1-cat snow depth
642      !!                zai : 1-cat ice concentration
643      !!
644      !! ** Output    : jpl-cat
645      !!
646      !!  (Example of application: BDY forcings when input are cell averaged) 
647      !!
648      !!-------------------------------------------------------------------
649      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
650      !!                    2014    (C. Rousset)        Rewriting
651      !!-------------------------------------------------------------------
652      !! Local variables
653      INTEGER  :: ji, jk, jl             ! dummy loop indices
654      INTEGER  :: ijpij, i_fill, jl0 
655      REAL(wp) :: zarg, zV, zconv, zdh, zdv
656      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
657      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
658      INTEGER , POINTER, DIMENSION(:)         ::   itest
659 
660      Call wrk_alloc( 4, itest )
661      !--------------------------------------------------------------------
662      ! initialisation of variables
663      !--------------------------------------------------------------------
664      ijpij = SIZE(zhti,1)
665      zht_i(1:ijpij,1:jpl) = 0._wp
666      zht_s(1:ijpij,1:jpl) = 0._wp
667      za_i (1:ijpij,1:jpl) = 0._wp
668
669      ! ----------------------------------------
670      ! distribution over the jpl ice categories
671      ! ----------------------------------------
672      DO ji = 1, ijpij
673         
674         IF( zhti(ji) > 0._wp ) THEN
675
676            ! find which category (jl0) the input ice thickness falls into
677            jl0 = jpl
678            DO jl = 1, jpl
679               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
680                  jl0 = jl
681                  CYCLE
682               ENDIF
683            END DO
684
685            ! initialisation of tests
686            itest(:)  = 0
687         
688            i_fill = jpl + 1                                             !====================================
689            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
690               ! iteration                                               !====================================
691               i_fill = i_fill - 1
692               
693               ! initialisation of ice variables for each try
694               zht_i(ji,1:jpl) = 0._wp
695               za_i (ji,1:jpl) = 0._wp
696               itest(:)        = 0           
697               
698               ! *** case very thin ice: fill only category 1
699               IF ( i_fill == 1 ) THEN
700                  zht_i(ji,1) = zhti(ji)
701                  za_i (ji,1) = zai (ji)
702                 
703               ! *** case ice is thicker: fill categories >1
704               ELSE
705
706                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean
707                  DO jl = 1, i_fill - 1
708                     zht_i(ji,jl) = hi_mean(jl)
709                  END DO
710                 
711                  ! Concentrations in the (i_fill-1) categories
712                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
713                  DO jl = 1, i_fill - 1
714                     IF ( jl /= jl0 ) THEN
715                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
716                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
717                     ENDIF
718                  END DO
719                 
720                  ! Concentration in the last (i_fill) category
721                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
722                 
723                  ! Ice thickness in the last (i_fill) category
724                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
725                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
726                 
727                  ! clem: correction if concentration of upper cat is greater than lower cat
728                  !       (it should be a gaussian around jl0 but sometimes it is not)
729                  IF ( jl0 /= jpl ) THEN
730                     DO jl = jpl, jl0+1, -1
731                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
732                           zdv = zht_i(ji,jl) * za_i(ji,jl)
733                           zht_i(ji,jl    ) = 0._wp
734                           za_i (ji,jl    ) = 0._wp
735                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
736                        END IF
737                     ENDDO
738                  ENDIF
739               
740               ENDIF ! case ice is thick or thin
741               
742               !---------------------
743               ! Compatibility tests
744               !---------------------
745               ! Test 1: area conservation
746               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
747               IF ( zconv < epsi06 ) itest(1) = 1
748           
749               ! Test 2: volume conservation
750               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
751               IF ( zconv < epsi06 ) itest(2) = 1
752               
753               ! Test 3: thickness of the last category is in-bounds ?
754               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
755               
756               ! Test 4: positivity of ice concentrations
757               itest(4) = 1
758               DO jl = 1, i_fill
759                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
760               END DO
761               !                                         !============================
762            END DO                                       ! end iteration on categories
763               !                                         !============================
764         ENDIF ! if zhti > 0
765      END DO ! i loop
766     
767      ! ------------------------------------------------
768      ! Adding Snow in each category where za_i is not 0
769      ! ------------------------------------------------
770      DO jl = 1, jpl
771         DO ji = 1, ijpij
772            IF( za_i(ji,jl) > 0._wp ) THEN
773               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
774               ! In case snow load is in excess that would lead to transformation from snow to ice
775               ! Then, transfer the snow excess into the ice (different from limthd_dh)
776               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
777               ! recompute ht_i, ht_s avoiding out of bounds values
778               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
779               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
780            ENDIF
781         ENDDO
782      ENDDO
783
784      CALL wrk_dealloc( 4, itest )
785      !
786    END SUBROUTINE lim_var_itd
787
788
789#else
790   !!----------------------------------------------------------------------
791   !!   Default option         Dummy module          NO  LIM3 sea-ice model
792   !!----------------------------------------------------------------------
793CONTAINS
794   SUBROUTINE lim_var_agg          ! Empty routines
795   END SUBROUTINE lim_var_agg
796   SUBROUTINE lim_var_glo2eqv      ! Empty routines
797   END SUBROUTINE lim_var_glo2eqv
798   SUBROUTINE lim_var_eqv2glo      ! Empty routines
799   END SUBROUTINE lim_var_eqv2glo
800   SUBROUTINE lim_var_salprof      ! Empty routines
801   END SUBROUTINE lim_var_salprof
802   SUBROUTINE lim_var_bv           ! Emtpy routines
803   END SUBROUTINE lim_var_bv
804   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
805   END SUBROUTINE lim_var_salprof1d
806   SUBROUTINE lim_var_zapsmall
807   END SUBROUTINE lim_var_zapsmall
808   SUBROUTINE lim_var_itd
809   END SUBROUTINE lim_var_itd
810#endif
811
812   !!======================================================================
813END MODULE limvar
Note: See TracBrowser for help on using the repository browser.