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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 5179

Last change on this file since 5179 was 5179, checked in by clem, 9 years ago

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