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

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

Commit a first set of modifications
  • Property svn:keywords set to Id
File size: 34.4 KB
Line 
1MODULE limvar
2   !!======================================================================
3   !!                       ***  MODULE limvar ***
4   !!                 Different sets of ice model variables
5   !!                   how to switch from one to another
6   !!
7   !!                 There are three sets of variables
8   !!                 VGLO : global variables of the model
9   !!                        - v_i (jpi,jpj,jpl)
10   !!                        - v_s (jpi,jpj,jpl)
11   !!                        - a_i (jpi,jpj,jpl)
12   !!                        - t_s (jpi,jpj,jpl)
13   !!                        - e_i (jpi,jpj,nlay_i,jpl)
14   !!                        - smv_i(jpi,jpj,jpl)
15   !!                        - oa_i (jpi,jpj,jpl)
16   !!                 VEQV : equivalent variables sometimes used in the model
17   !!                        - ht_i(jpi,jpj,jpl)
18   !!                        - ht_s(jpi,jpj,jpl)
19   !!                        - t_i (jpi,jpj,nlay_i,jpl)
20   !!                        ...
21   !!                 VAGG : aggregate variables, averaged/summed over all
22   !!                        thickness categories
23   !!                        - vt_i(jpi,jpj)
24   !!                        - vt_s(jpi,jpj)
25   !!                        - at_i(jpi,jpj)
26   !!                        - et_s(jpi,jpj)  !total snow heat content
27   !!                        - et_i(jpi,jpj)  !total ice thermal content
28   !!                        - smt_i(jpi,jpj) !mean ice salinity
29   !!                        - 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 ! MV SIMIP 2016
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               ! SIMIP snow melt diagnostic
603               diag_dms_mel_1d(ji) = diag_dms_mel(ji,jj) + ( v_s(ji,jj,jl) - zvs  ) * rhosn * r1_rdtice
604
605            END DO
606         END DO
607      END DO 
608
609      ! to be sure that at_i is the sum of a_i(jl)
610      at_i (:,:) = 0._wp
611      DO jl = 1, jpl
612         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
613      END DO
614
615      ! open water = 1 if at_i=0
616      DO jj = 1, jpj
617         DO ji = 1, jpi
618            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
619            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
620         END DO
621      END DO
622
623      !
624   END SUBROUTINE lim_var_zapsmall
625
626   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
627      !!------------------------------------------------------------------
628      !!                ***  ROUTINE lim_var_itd   ***
629      !!
630      !! ** Purpose :  converting 1-cat ice to multiple ice categories
631      !!
632      !!                  ice thickness distribution follows a gaussian law
633      !!               around the concentration of the most likely ice thickness
634      !!                           (similar as limistate.F90)
635      !!
636      !! ** Method:   Iterative procedure
637      !!               
638      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
639      !!
640      !!               2) Check whether the distribution conserves area and volume, positivity and
641      !!                  category boundaries
642      !!             
643      !!               3) If not (input ice is too thin), the last category is empty and
644      !!                  the number of categories is reduced (jpl-1)
645      !!
646      !!               4) Iterate until ok (SUM(itest(:) = 4)
647      !!
648      !! ** Arguments : zhti: 1-cat ice thickness
649      !!                zhts: 1-cat snow depth
650      !!                zai : 1-cat ice concentration
651      !!
652      !! ** Output    : jpl-cat
653      !!
654      !!  (Example of application: BDY forcings when input are cell averaged) 
655      !!
656      !!-------------------------------------------------------------------
657      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
658      !!                    2014    (C. Rousset)        Rewriting
659      !!-------------------------------------------------------------------
660      !! Local variables
661      INTEGER  :: ji, jk, jl             ! dummy loop indices
662      INTEGER  :: ijpij, i_fill, jl0 
663      REAL(wp) :: zarg, zV, zconv, zdh
664      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
665      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
666      INTEGER , POINTER, DIMENSION(:)         ::   itest
667 
668      CALL wrk_alloc( 4, itest )
669      !--------------------------------------------------------------------
670      ! initialisation of variables
671      !--------------------------------------------------------------------
672      ijpij = SIZE(zhti,1)
673      zht_i(1:ijpij,1:jpl) = 0._wp
674      zht_s(1:ijpij,1:jpl) = 0._wp
675      za_i (1:ijpij,1:jpl) = 0._wp
676
677      ! ----------------------------------------
678      ! distribution over the jpl ice categories
679      ! ----------------------------------------
680      DO ji = 1, ijpij
681         
682         IF( zhti(ji) > 0._wp ) THEN
683
684         ! initialisation of tests
685         itest(:)  = 0
686         
687         i_fill = jpl + 1                                             !====================================
688         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
689            ! iteration                                               !====================================
690            i_fill = i_fill - 1
691           
692            ! initialisation of ice variables for each try
693            zht_i(ji,1:jpl) = 0._wp
694            za_i (ji,1:jpl) = 0._wp
695            itest(:)        = 0     
696     
697            ! *** case very thin ice: fill only category 1
698            IF ( i_fill == 1 ) THEN
699               zht_i(ji,1) = zhti(ji)
700               za_i (ji,1) = zai (ji)
701
702            ! *** case ice is thicker: fill categories >1
703            ELSE
704
705               ! Fill ice thicknesses except the last one (i_fill) by hmean
706               DO jl = 1, i_fill - 1
707                  zht_i(ji,jl) = hi_mean(jl)
708               END DO
709               
710               ! find which category (jl0) the input ice thickness falls into
711               jl0 = i_fill
712               DO jl = 1, i_fill
713                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
714                     jl0 = jl
715           CYCLE
716                  ENDIF
717               END DO
718               
719               ! Concentrations in the (i_fill-1) categories
720               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
721               DO jl = 1, i_fill - 1
722                  IF ( jl == jl0 ) CYCLE
723                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
724                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
725               END DO
726               
727               ! Concentration in the last (i_fill) category
728               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
729               
730               ! Ice thickness in the last (i_fill) category
731               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
732               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
733               
734            ENDIF ! case ice is thick or thin
735           
736            !---------------------
737            ! Compatibility tests
738            !---------------------
739            ! Test 1: area conservation
740            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
741            IF ( zconv < epsi06 ) itest(1) = 1
742           
743            ! Test 2: volume conservation
744            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
745            IF ( zconv < epsi06 ) itest(2) = 1
746           
747            ! Test 3: thickness of the last category is in-bounds ?
748            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
749           
750            ! Test 4: positivity of ice concentrations
751            itest(4) = 1
752            DO jl = 1, i_fill
753               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
754            END DO           
755                                                           !============================
756         END DO                                            ! end iteration on categories
757                                                           !============================
758         ENDIF ! if zhti > 0
759      END DO ! i loop
760
761      ! ------------------------------------------------
762      ! Adding Snow in each category where za_i is not 0
763      ! ------------------------------------------------
764      DO jl = 1, jpl
765         DO ji = 1, ijpij
766            IF( za_i(ji,jl) > 0._wp ) THEN
767               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
768               ! In case snow load is in excess that would lead to transformation from snow to ice
769               ! Then, transfer the snow excess into the ice (different from limthd_dh)
770               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
771               ! recompute ht_i, ht_s avoiding out of bounds values
772               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
773               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
774            ENDIF
775         ENDDO
776      ENDDO
777
778      CALL wrk_dealloc( 4, itest )
779      !
780    END SUBROUTINE lim_var_itd
781
782
783#else
784   !!----------------------------------------------------------------------
785   !!   Default option         Dummy module          NO  LIM3 sea-ice model
786   !!----------------------------------------------------------------------
787CONTAINS
788   SUBROUTINE lim_var_agg          ! Empty routines
789   END SUBROUTINE lim_var_agg
790   SUBROUTINE lim_var_glo2eqv      ! Empty routines
791   END SUBROUTINE lim_var_glo2eqv
792   SUBROUTINE lim_var_eqv2glo      ! Empty routines
793   END SUBROUTINE lim_var_eqv2glo
794   SUBROUTINE lim_var_salprof      ! Empty routines
795   END SUBROUTINE lim_var_salprof
796   SUBROUTINE lim_var_bv           ! Emtpy routines
797   END SUBROUTINE lim_var_bv
798   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
799   END SUBROUTINE lim_var_salprof1d
800   SUBROUTINE lim_var_zapsmall
801   END SUBROUTINE lim_var_zapsmall
802   SUBROUTINE lim_var_itd
803   END SUBROUTINE lim_var_itd
804#endif
805
806   !!======================================================================
807END MODULE limvar
Note: See TracBrowser for help on using the repository browser.