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

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

All wrk_alloc removed

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