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

source: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 8152

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

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