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

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

LIM3 change all namelist names to fit with NEMO convention
  • Property svn:keywords set to Id
File size: 32.6 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) + rtt              ! Ice layer melt temperature
196                  !
197                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
198                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i / rhoic - lfus
199                  zccc       =  lfus * (ztmelts-rtt)
200                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
201                  t_i(ji,jj,jk,jl) = rtt + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )
202                  t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rtt < t_i < rtt
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) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 )
222                  t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rtt < t_i < rtt
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) ::   dummy_fac0, dummy_fac1, dummy_fac, zsal
283      REAL(wp) ::   zswi0, zswi01, zswibal, 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         dummy_fac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
314         dummy_fac1 = 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 zswibal = 1
325                  ! this is to force a constant salinity profile in the Baltic Sea
326                  zswibal = 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) * dummy_fac0 + dummy_fac1 )
328                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal )
329               END DO
330            END DO
331         END DO
332
333         dummy_fac = 1._wp / REAL( nlay_i )                   ! 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) * dummy_fac
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 ) / REAL(nlay_i,wp)
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                     &                      / (  REAL(nlay_i,wp) * 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) - rtt) + epsi10 ) )  )
418                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   &
419                     &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp)
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) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars
442      REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, 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         dummy_fac0 = 1._wp / ( zsi0 - zsi1 )
469         dummy_fac1 = zsi1 / ( zsi1 - zsi0 )
470         dummy_fac2 = 1._wp / REAL(nlay_i,wp)
471
472         DO jk = 1, nlay_i
473            DO ji = kideb, kiut
474               ii =  MOD( npb(ji) - 1 , jpi ) + 1
475               ij =     ( npb(ji) - 1 ) / jpi + 1
476               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
477               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
478               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
479               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
480               ! if 2.sm_i GE sss_m then zswibal = 1
481               ! this is to force a constant salinity profile in the Baltic Sea
482               zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
483               !
484               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zswibal )
485               !
486               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2
487               ! weighting the profile
488               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
489            END DO
490         END DO
491
492      ENDIF 
493
494      !-------------------------------------------------------
495      ! Vertically varying salinity profile, constant in time
496      !-------------------------------------------------------
497
498      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
499         !
500         sm_i_1d(:) = 2.30_wp
501         !
502         DO jk = 1, nlay_i
503            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
504            zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
505            DO ji = kideb, kiut
506               s_i_1d(ji,jk) = zsal
507            END DO
508         END DO
509         !
510      ENDIF
511      !
512      CALL wrk_dealloc( jpij, z_slope_s )
513      !
514   END SUBROUTINE lim_var_salprof1d
515
516   SUBROUTINE lim_var_zapsmall
517      !!-------------------------------------------------------------------
518      !!                   ***  ROUTINE lim_var_zapsmall ***
519      !!
520      !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes
521      !!
522      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
523      !!-------------------------------------------------------------------
524      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
525
526      REAL(wp) ::   zsal, zvi, zvs, zei, zes
527      !!-------------------------------------------------------------------
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                  zei              = e_i(ji,jj,jk,jl)
539                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
540                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rtt * ( 1._wp - rswitch )
541                  ! update exchanges with ocean
542                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
543               END DO
544            END DO
545         END DO
546
547         DO jj = 1 , jpj
548            DO ji = 1 , jpi
549               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
550               
551               zsal = smv_i(ji,jj,  jl)
552               zvi  = v_i  (ji,jj,  jl)
553               zvs  = v_s  (ji,jj,  jl)
554               zes  = e_s  (ji,jj,1,jl)
555               !-----------------------------------------------------------------
556               ! Zap snow energy
557               !-----------------------------------------------------------------
558               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rtt * ( 1._wp - rswitch )
559               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
560
561               !-----------------------------------------------------------------
562               ! zap ice and snow volume, add water and salt to ocean
563               !-----------------------------------------------------------------
564               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
565               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
566               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
567               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
568               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
569               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
570               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
571
572               ! ice salinity must stay in bounds
573               IF(  nn_icesal == 2  ) THEN
574                  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) )
575               ENDIF
576               ! update exchanges with ocean
577               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
578               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
579               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
580               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
581            END DO
582         END DO
583      END DO ! jl
584
585      ! to be sure that at_i is the sum of a_i(jl)
586      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
587      !
588   END SUBROUTINE lim_var_zapsmall
589
590   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
591      !!------------------------------------------------------------------
592      !!                ***  ROUTINE lim_var_itd   ***
593      !!
594      !! ** Purpose :  converting 1-cat ice to multiple ice categories
595      !!
596      !!                  ice thickness distribution follows a gaussian law
597      !!               around the concentration of the most likely ice thickness
598      !!                           (similar as limistate.F90)
599      !!
600      !! ** Method:   Iterative procedure
601      !!               
602      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
603      !!
604      !!               2) Check whether the distribution conserves area and volume, positivity and
605      !!                  category boundaries
606      !!             
607      !!               3) If not (input ice is too thin), the last category is empty and
608      !!                  the number of categories is reduced (jpl-1)
609      !!
610      !!               4) Iterate until ok (SUM(itest(:) = 4)
611      !!
612      !! ** Arguments : zhti: 1-cat ice thickness
613      !!                zhts: 1-cat snow depth
614      !!                zai : 1-cat ice concentration
615      !!
616      !! ** Output    : jpl-cat
617      !!
618      !!  (Example of application: BDY forcings when input are cell averaged) 
619      !!
620      !!-------------------------------------------------------------------
621      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
622      !!                    2014    (C. Rousset)        Rewriting
623      !!-------------------------------------------------------------------
624      !! Local variables
625      INTEGER  :: ji, jk, jl             ! dummy loop indices
626      INTEGER  :: ijpij, i_fill, jl0 
627      REAL(wp) :: zarg, zV, zconv, zdh
628      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
629      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
630      INTEGER , POINTER, DIMENSION(:)         ::   itest
631 
632      CALL wrk_alloc( 4, itest )
633      !--------------------------------------------------------------------
634      ! initialisation of variables
635      !--------------------------------------------------------------------
636      ijpij = SIZE(zhti,1)
637      zht_i(1:ijpij,1:jpl) = 0._wp
638      zht_s(1:ijpij,1:jpl) = 0._wp
639      za_i (1:ijpij,1:jpl) = 0._wp
640
641      ! ----------------------------------------
642      ! distribution over the jpl ice categories
643      ! ----------------------------------------
644      DO ji = 1, ijpij
645         
646         IF( zhti(ji) > 0._wp ) THEN
647
648         ! initialisation of tests
649         itest(:)  = 0
650         
651         i_fill = jpl + 1                                             !====================================
652         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
653            ! iteration                                               !====================================
654            i_fill = i_fill - 1
655           
656            ! initialisation of ice variables for each try
657            zht_i(ji,1:jpl) = 0._wp
658            za_i (ji,1:jpl) = 0._wp
659           
660            ! *** case very thin ice: fill only category 1
661            IF ( i_fill == 1 ) THEN
662               zht_i(ji,1) = zhti(ji)
663               za_i (ji,1) = zai (ji)
664
665            ! *** case ice is thicker: fill categories >1
666            ELSE
667
668               ! Fill ice thicknesses except the last one (i_fill) by hmean
669               DO jl = 1, i_fill - 1
670                  zht_i(ji,jl) = hi_mean(jl)
671               END DO
672               
673               ! find which category (jl0) the input ice thickness falls into
674               jl0 = i_fill
675               DO jl = 1, i_fill
676                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
677                     jl0 = jl
678           CYCLE
679                  ENDIF
680               END DO
681               
682               ! Concentrations in the (i_fill-1) categories
683               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
684               DO jl = 1, i_fill - 1
685                  IF ( jl == jl0 ) CYCLE
686                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
687                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
688               END DO
689               
690               ! Concentration in the last (i_fill) category
691               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
692               
693               ! Ice thickness in the last (i_fill) category
694               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
695               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
696               
697            ENDIF ! case ice is thick or thin
698           
699            !---------------------
700            ! Compatibility tests
701            !---------------------
702            ! Test 1: area conservation
703            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
704            IF ( zconv < epsi06 ) itest(1) = 1
705           
706            ! Test 2: volume conservation
707            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
708            IF ( zconv < epsi06 ) itest(2) = 1
709           
710            ! Test 3: thickness of the last category is in-bounds ?
711            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
712           
713            ! Test 4: positivity of ice concentrations
714            itest(4) = 1
715            DO jl = 1, i_fill
716               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
717            END DO           
718                                                           !============================
719         END DO                                            ! end iteration on categories
720                                                           !============================
721         ENDIF ! if zhti > 0
722      END DO ! i loop
723
724      ! ------------------------------------------------
725      ! Adding Snow in each category where za_i is not 0
726      ! ------------------------------------------------
727      DO jl = 1, jpl
728         DO ji = 1, ijpij
729            IF( za_i(ji,jl) > 0._wp ) THEN
730               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
731               ! In case snow load is in excess that would lead to transformation from snow to ice
732               ! Then, transfer the snow excess into the ice (different from limthd_dh)
733               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
734               ! recompute ht_i, ht_s avoiding out of bounds values
735               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
736               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
737            ENDIF
738         ENDDO
739      ENDDO
740
741      CALL wrk_dealloc( 4, itest )
742      !
743    END SUBROUTINE lim_var_itd
744
745
746#else
747   !!----------------------------------------------------------------------
748   !!   Default option         Dummy module          NO  LIM3 sea-ice model
749   !!----------------------------------------------------------------------
750CONTAINS
751   SUBROUTINE lim_var_agg          ! Empty routines
752   END SUBROUTINE lim_var_agg
753   SUBROUTINE lim_var_glo2eqv      ! Empty routines
754   END SUBROUTINE lim_var_glo2eqv
755   SUBROUTINE lim_var_eqv2glo      ! Empty routines
756   END SUBROUTINE lim_var_eqv2glo
757   SUBROUTINE lim_var_salprof      ! Empty routines
758   END SUBROUTINE lim_var_salprof
759   SUBROUTINE lim_var_bv           ! Emtpy routines
760   END SUBROUTINE lim_var_bv
761   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
762   END SUBROUTINE lim_var_salprof1d
763   SUBROUTINE lim_var_zapsmall
764   END SUBROUTINE lim_var_zapsmall
765   SUBROUTINE lim_var_itd
766   END SUBROUTINE lim_var_itd
767#endif
768
769   !!======================================================================
770END MODULE limvar
Note: See TracBrowser for help on using the repository browser.