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

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 7560

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

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