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

Last change on this file since 7293 was 7293, checked in by vancop, 7 years ago

Commit infrastructure for melt pond impletation

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