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

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

major LIM3 cleaning + monocat capabilities + NEMO namelist-consistency; sette to follow
  • Property svn:keywords set to Id
File size: 33.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   !!                        - 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) - epsi10 ) ) 
127                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch   ! ice salinity
128                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
129                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * 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 = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !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) , epsi10 ) * rswitch
165               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
166               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * 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 = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !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) , epsi10 ) * 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   = 1.0 - 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 = 1._wp - 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      tm_i(:,:) = 0._wp
231      DO jl = 1, jpl
232         DO jk = 1, nlay_i
233            DO jj = 1, jpj
234               DO ji = 1, jpi
235                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
236                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
237                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
238               END DO
239            END DO
240         END DO
241      END DO
242      !
243   END SUBROUTINE lim_var_glo2eqv
244
245
246   SUBROUTINE lim_var_eqv2glo
247      !!------------------------------------------------------------------
248      !!                ***  ROUTINE lim_var_eqv2glo ***
249      !!
250      !! ** Purpose :   computes global variables as function of equivalent variables
251      !!                i.e. it turns VEQV into VGLO
252      !! ** Method  :
253      !!
254      !! ** History :  (01-2006) Martin Vancoppenolle, UCL-ASTR
255      !!------------------------------------------------------------------
256      !
257      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
258      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
259      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
260      oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:)
261      !
262   END SUBROUTINE lim_var_eqv2glo
263
264
265   SUBROUTINE lim_var_salprof
266      !!------------------------------------------------------------------
267      !!                ***  ROUTINE lim_var_salprof ***
268      !!
269      !! ** Purpose :   computes salinity profile in function of bulk salinity     
270      !!
271      !! ** Method  : If bulk salinity greater than zsi1,
272      !!              the profile is assumed to be constant (S_inf)
273      !!              If bulk salinity lower than zsi0,
274      !!              the profile is linear with 0 at the surface (S_zero)
275      !!              If it is between zsi0 and zsi1, it is a
276      !!              alpha-weighted linear combination of s_inf and s_zero
277      !!
278      !! ** References : Vancoppenolle et al., 2007
279      !!------------------------------------------------------------------
280      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
281      REAL(wp) ::   zfac0, zfac1, zsal
282      REAL(wp) ::   zswi0, zswi01, zargtemp , zs_zero   
283      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha
284      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
285      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
286      !!------------------------------------------------------------------
287
288      CALL wrk_alloc( jpi, jpj, jpl, z_slope_s, zalpha )
289
290      !---------------------------------------
291      ! Vertically constant, constant in time
292      !---------------------------------------
293      IF(  nn_icesal == 1  )   s_i(:,:,:,:) = rn_icesal
294
295      !-----------------------------------
296      ! Salinity profile, varying in time
297      !-----------------------------------
298      IF(  nn_icesal == 2  ) THEN
299         !
300         DO jk = 1, nlay_i
301            s_i(:,:,jk,:)  = sm_i(:,:,:)
302         END DO
303         !
304         DO jl = 1, jpl                               ! Slope of the linear profile
305            DO jj = 1, jpj
306               DO ji = 1, jpi
307                  z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) )
308               END DO
309            END DO
310         END DO
311         !
312         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
313         zfac1 = zsi1  / ( zsi1 - zsi0 )
314         !
315         zalpha(:,:,:) = 0._wp
316         DO jl = 1, jpl
317            DO jj = 1, jpj
318               DO ji = 1, jpi
319                  ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
320                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , zsi0 - sm_i(ji,jj,jl) ) ) 
321                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
322                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) ) 
323                  ! If 2.sm_i GE sss_m then rswitch = 1
324                  ! this is to force a constant salinity profile in the Baltic Sea
325                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
326                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
327                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
328               END DO
329            END DO
330         END DO
331
332         ! Computation of the profile
333         DO jl = 1, jpl
334            DO jk = 1, nlay_i
335               DO jj = 1, jpj
336                  DO ji = 1, jpi
337                     !                                      ! linear profile with 0 at the surface
338                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
339                     !                                      ! weighting the profile
340                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
341                  END DO ! ji
342               END DO ! jj
343            END DO ! jk
344         END DO ! jl
345         !
346      ENDIF ! nn_icesal
347
348      !-------------------------------------------------------
349      ! Vertically varying salinity profile, constant in time
350      !-------------------------------------------------------
351
352      IF(  nn_icesal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
353         !
354         sm_i(:,:,:) = 2.30_wp
355         !
356         DO jl = 1, jpl
357            DO jk = 1, nlay_i
358               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
359               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
360               s_i(:,:,jk,jl) =  zsal
361            END DO
362         END DO
363         !
364      ENDIF ! nn_icesal
365      !
366      CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
367      !
368   END SUBROUTINE lim_var_salprof
369
370
371   SUBROUTINE lim_var_icetm
372      !!------------------------------------------------------------------
373      !!                ***  ROUTINE lim_var_icetm ***
374      !!
375      !! ** Purpose :   computes mean sea ice temperature
376      !!------------------------------------------------------------------
377      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
378      !!------------------------------------------------------------------
379
380      ! Mean sea ice temperature
381      tm_i(:,:) = 0._wp
382      DO jl = 1, jpl
383         DO jk = 1, nlay_i
384            DO jj = 1, jpj
385               DO ji = 1, jpi
386                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
387                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
388                     &                      * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 )
389               END DO
390            END DO
391         END DO
392      END DO
393
394   END SUBROUTINE lim_var_icetm
395
396
397   SUBROUTINE lim_var_bv
398      !!------------------------------------------------------------------
399      !!                ***  ROUTINE lim_var_bv ***
400      !!
401      !! ** Purpose :   computes mean brine volume (%) in sea ice
402      !!
403      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
404      !!
405      !! References : Vancoppenolle et al., JGR, 2007
406      !!------------------------------------------------------------------
407      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
408      REAL(wp) ::   zbvi             ! local scalars
409      !!------------------------------------------------------------------
410      !
411      bv_i(:,:) = 0._wp
412      DO jl = 1, jpl
413         DO jk = 1, nlay_i
414            DO jj = 1, jpj
415               DO ji = 1, jpi
416                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
417                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   &
418                     &                   * v_i(ji,jj,jl) * r1_nlay_i
419                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
420                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi10 )
421               END DO
422            END DO
423         END DO
424      END DO
425      !
426   END SUBROUTINE lim_var_bv
427
428
429   SUBROUTINE lim_var_salprof1d( kideb, kiut )
430      !!-------------------------------------------------------------------
431      !!                  ***  ROUTINE lim_thd_salprof1d  ***
432      !!
433      !! ** Purpose :   1d computation of the sea ice salinity profile
434      !!                Works with 1d vectors and is used by thermodynamic modules
435      !!-------------------------------------------------------------------
436      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
437      !
438      INTEGER  ::   ji, jk    ! dummy loop indices
439      INTEGER  ::   ii, ij    ! local integers
440      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars
441      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      -
442      !
443      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
444      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
445      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
446      !!---------------------------------------------------------------------
447
448      CALL wrk_alloc( jpij, z_slope_s )
449
450      !---------------------------------------
451      ! Vertically constant, constant in time
452      !---------------------------------------
453      IF( nn_icesal == 1 )   s_i_1d(:,:) = rn_icesal
454
455      !------------------------------------------------------
456      ! Vertically varying salinity profile, varying in time
457      !------------------------------------------------------
458
459      IF(  nn_icesal == 2  ) THEN
460         !
461         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
462            z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) )
463         END DO
464
465         ! Weighting factor between zs_zero and zs_inf
466         !---------------------------------------------
467         zfac0 = 1._wp / ( zsi0 - zsi1 )
468         zfac1 = zsi1 / ( zsi1 - zsi0 )
469         DO jk = 1, nlay_i
470            DO ji = kideb, kiut
471               ii =  MOD( npb(ji) - 1 , jpi ) + 1
472               ij =     ( npb(ji) - 1 ) / jpi + 1
473               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
474               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
475               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
476               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
477               ! if 2.sm_i GE sss_m then rswitch = 1
478               ! this is to force a constant salinity profile in the Baltic Sea
479               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
480               !
481               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch )
482               !
483               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
484               ! weighting the profile
485               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
486            END DO
487         END DO
488
489      ENDIF 
490
491      !-------------------------------------------------------
492      ! Vertically varying salinity profile, constant in time
493      !-------------------------------------------------------
494
495      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
496         !
497         sm_i_1d(:) = 2.30_wp
498         !
499         DO jk = 1, nlay_i
500            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
501            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
502            DO ji = kideb, kiut
503               s_i_1d(ji,jk) = zsal
504            END DO
505         END DO
506         !
507      ENDIF
508      !
509      CALL wrk_dealloc( jpij, z_slope_s )
510      !
511   END SUBROUTINE lim_var_salprof1d
512
513   SUBROUTINE lim_var_zapsmall
514      !!-------------------------------------------------------------------
515      !!                   ***  ROUTINE lim_var_zapsmall ***
516      !!
517      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
518      !!
519      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
520      !!-------------------------------------------------------------------
521      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
522      REAL(wp) ::   zsal, zvi, zvs, zei, zes
523      !!-------------------------------------------------------------------
524      at_i (:,:) = 0._wp
525      DO jl = 1, jpl
526         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
527      END DO
528
529      DO jl = 1, jpl
530
531         !-----------------------------------------------------------------
532         ! Zap ice energy and use ocean heat to melt ice
533         !-----------------------------------------------------------------
534         DO jk = 1, nlay_i
535            DO jj = 1 , jpj
536               DO ji = 1 , jpi
537                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
538                  rswitch          = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
539                  zei              = e_i(ji,jj,jk,jl)
540                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
541                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
542                  ! update exchanges with ocean
543                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
544               END DO
545            END DO
546         END DO
547
548         DO jj = 1 , jpj
549            DO ji = 1 , jpi
550               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
551               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch
552               
553               zsal = smv_i(ji,jj,  jl)
554               zvi  = v_i  (ji,jj,  jl)
555               zvs  = v_s  (ji,jj,  jl)
556               zes  = e_s  (ji,jj,1,jl)
557               !-----------------------------------------------------------------
558               ! Zap snow energy
559               !-----------------------------------------------------------------
560               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
561               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
562
563               !-----------------------------------------------------------------
564               ! zap ice and snow volume, add water and salt to ocean
565               !-----------------------------------------------------------------
566               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
567               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
568               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
569               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
570               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
571               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
572               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
573
574               ! ice salinity must stay in bounds
575               IF(  nn_icesal == 2  ) THEN
576                  smv_i(ji,jj,jl) = MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) )
577               ENDIF
578               ! update exchanges with ocean
579               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
580               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
581               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
582               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
583            END DO
584         END DO
585      END DO 
586
587      ! to be sure that at_i is the sum of a_i(jl)
588      at_i (:,:) = 0._wp
589      DO jl = 1, jpl
590         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
591      END DO
592
593      ! open water = 1 if at_i=0
594      DO jj = 1, jpj
595         DO ji = 1, jpi
596            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
597            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
598         END DO
599      END DO
600
601      !
602   END SUBROUTINE lim_var_zapsmall
603
604   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
605      !!------------------------------------------------------------------
606      !!                ***  ROUTINE lim_var_itd   ***
607      !!
608      !! ** Purpose :  converting 1-cat ice to multiple ice categories
609      !!
610      !!                  ice thickness distribution follows a gaussian law
611      !!               around the concentration of the most likely ice thickness
612      !!                           (similar as limistate.F90)
613      !!
614      !! ** Method:   Iterative procedure
615      !!               
616      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
617      !!
618      !!               2) Check whether the distribution conserves area and volume, positivity and
619      !!                  category boundaries
620      !!             
621      !!               3) If not (input ice is too thin), the last category is empty and
622      !!                  the number of categories is reduced (jpl-1)
623      !!
624      !!               4) Iterate until ok (SUM(itest(:) = 4)
625      !!
626      !! ** Arguments : zhti: 1-cat ice thickness
627      !!                zhts: 1-cat snow depth
628      !!                zai : 1-cat ice concentration
629      !!
630      !! ** Output    : jpl-cat
631      !!
632      !!  (Example of application: BDY forcings when input are cell averaged) 
633      !!
634      !!-------------------------------------------------------------------
635      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
636      !!                    2014    (C. Rousset)        Rewriting
637      !!-------------------------------------------------------------------
638      !! Local variables
639      INTEGER  :: ji, jk, jl             ! dummy loop indices
640      INTEGER  :: ijpij, i_fill, jl0 
641      REAL(wp) :: zarg, zV, zconv, zdh
642      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
643      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
644      INTEGER , POINTER, DIMENSION(:)         ::   itest
645 
646      CALL wrk_alloc( 4, itest )
647      !--------------------------------------------------------------------
648      ! initialisation of variables
649      !--------------------------------------------------------------------
650      ijpij = SIZE(zhti,1)
651      zht_i(1:ijpij,1:jpl) = 0._wp
652      zht_s(1:ijpij,1:jpl) = 0._wp
653      za_i (1:ijpij,1:jpl) = 0._wp
654
655      ! ----------------------------------------
656      ! distribution over the jpl ice categories
657      ! ----------------------------------------
658      DO ji = 1, ijpij
659         
660         IF( zhti(ji) > 0._wp ) THEN
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           
674            ! *** case very thin ice: fill only category 1
675            IF ( i_fill == 1 ) THEN
676               zht_i(ji,1) = zhti(ji)
677               za_i (ji,1) = zai (ji)
678
679            ! *** case ice is thicker: fill categories >1
680            ELSE
681
682               ! Fill ice thicknesses except the last one (i_fill) by hmean
683               DO jl = 1, i_fill - 1
684                  zht_i(ji,jl) = hi_mean(jl)
685               END DO
686               
687               ! find which category (jl0) the input ice thickness falls into
688               jl0 = i_fill
689               DO jl = 1, i_fill
690                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
691                     jl0 = jl
692           CYCLE
693                  ENDIF
694               END DO
695               
696               ! Concentrations in the (i_fill-1) categories
697               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
698               DO jl = 1, i_fill - 1
699                  IF ( jl == jl0 ) CYCLE
700                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
701                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
702               END DO
703               
704               ! Concentration in the last (i_fill) category
705               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
706               
707               ! Ice thickness in the last (i_fill) category
708               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
709               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
710               
711            ENDIF ! case ice is thick or thin
712           
713            !---------------------
714            ! Compatibility tests
715            !---------------------
716            ! Test 1: area conservation
717            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
718            IF ( zconv < epsi06 ) itest(1) = 1
719           
720            ! Test 2: volume conservation
721            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
722            IF ( zconv < epsi06 ) itest(2) = 1
723           
724            ! Test 3: thickness of the last category is in-bounds ?
725            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
726           
727            ! Test 4: positivity of ice concentrations
728            itest(4) = 1
729            DO jl = 1, i_fill
730               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
731            END DO           
732                                                           !============================
733         END DO                                            ! end iteration on categories
734                                                           !============================
735         ENDIF ! if zhti > 0
736      END DO ! i loop
737
738      ! ------------------------------------------------
739      ! Adding Snow in each category where za_i is not 0
740      ! ------------------------------------------------
741      DO jl = 1, jpl
742         DO ji = 1, ijpij
743            IF( za_i(ji,jl) > 0._wp ) THEN
744               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
745               ! In case snow load is in excess that would lead to transformation from snow to ice
746               ! Then, transfer the snow excess into the ice (different from limthd_dh)
747               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
748               ! recompute ht_i, ht_s avoiding out of bounds values
749               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
750               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
751            ENDIF
752         ENDDO
753      ENDDO
754
755      CALL wrk_dealloc( 4, itest )
756      !
757    END SUBROUTINE lim_var_itd
758
759
760#else
761   !!----------------------------------------------------------------------
762   !!   Default option         Dummy module          NO  LIM3 sea-ice model
763   !!----------------------------------------------------------------------
764CONTAINS
765   SUBROUTINE lim_var_agg          ! Empty routines
766   END SUBROUTINE lim_var_agg
767   SUBROUTINE lim_var_glo2eqv      ! Empty routines
768   END SUBROUTINE lim_var_glo2eqv
769   SUBROUTINE lim_var_eqv2glo      ! Empty routines
770   END SUBROUTINE lim_var_eqv2glo
771   SUBROUTINE lim_var_salprof      ! Empty routines
772   END SUBROUTINE lim_var_salprof
773   SUBROUTINE lim_var_bv           ! Emtpy routines
774   END SUBROUTINE lim_var_bv
775   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
776   END SUBROUTINE lim_var_salprof1d
777   SUBROUTINE lim_var_zapsmall
778   END SUBROUTINE lim_var_zapsmall
779   SUBROUTINE lim_var_itd
780   END SUBROUTINE lim_var_itd
781#endif
782
783   !!======================================================================
784END MODULE limvar
Note: See TracBrowser for help on using the repository browser.