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

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

LIM3: cosmetic changes to increase readability and performance
  • Property svn:keywords set to Id
File size: 32.4 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 * r1_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) ::   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) - 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) * 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.0 - 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 salt fluxes
519      !!
520      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
521      !!-------------------------------------------------------------------
522      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
523
524      REAL(wp) ::   zsal, zvi, zvs, zei, zes
525      !!-------------------------------------------------------------------
526
527      DO jl = 1, jpl
528
529         !-----------------------------------------------------------------
530         ! Zap ice energy and use ocean heat to melt ice
531         !-----------------------------------------------------------------
532         DO jk = 1, nlay_i
533            DO jj = 1 , jpj
534               DO ji = 1 , jpi
535                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
536                  zei              = e_i(ji,jj,jk,jl)
537                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
538                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rtt * ( 1._wp - rswitch )
539                  ! update exchanges with ocean
540                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
541               END DO
542            END DO
543         END DO
544
545         DO jj = 1 , jpj
546            DO ji = 1 , jpi
547               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
548               
549               zsal = smv_i(ji,jj,  jl)
550               zvi  = v_i  (ji,jj,  jl)
551               zvs  = v_s  (ji,jj,  jl)
552               zes  = e_s  (ji,jj,1,jl)
553               !-----------------------------------------------------------------
554               ! Zap snow energy
555               !-----------------------------------------------------------------
556               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rtt * ( 1._wp - rswitch )
557               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
558
559               !-----------------------------------------------------------------
560               ! zap ice and snow volume, add water and salt to ocean
561               !-----------------------------------------------------------------
562               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
563               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
564               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
565               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
566               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
567               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
568               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
569
570               ! ice salinity must stay in bounds
571               IF(  nn_icesal == 2  ) THEN
572                  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) )
573               ENDIF
574               ! update exchanges with ocean
575               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
576               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
577               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
578               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
579            END DO
580         END DO
581      END DO ! jl
582
583      ! to be sure that at_i is the sum of a_i(jl)
584      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
585      !
586   END SUBROUTINE lim_var_zapsmall
587
588   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
589      !!------------------------------------------------------------------
590      !!                ***  ROUTINE lim_var_itd   ***
591      !!
592      !! ** Purpose :  converting 1-cat ice to multiple ice categories
593      !!
594      !!                  ice thickness distribution follows a gaussian law
595      !!               around the concentration of the most likely ice thickness
596      !!                           (similar as limistate.F90)
597      !!
598      !! ** Method:   Iterative procedure
599      !!               
600      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
601      !!
602      !!               2) Check whether the distribution conserves area and volume, positivity and
603      !!                  category boundaries
604      !!             
605      !!               3) If not (input ice is too thin), the last category is empty and
606      !!                  the number of categories is reduced (jpl-1)
607      !!
608      !!               4) Iterate until ok (SUM(itest(:) = 4)
609      !!
610      !! ** Arguments : zhti: 1-cat ice thickness
611      !!                zhts: 1-cat snow depth
612      !!                zai : 1-cat ice concentration
613      !!
614      !! ** Output    : jpl-cat
615      !!
616      !!  (Example of application: BDY forcings when input are cell averaged) 
617      !!
618      !!-------------------------------------------------------------------
619      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
620      !!                    2014    (C. Rousset)        Rewriting
621      !!-------------------------------------------------------------------
622      !! Local variables
623      INTEGER  :: ji, jk, jl             ! dummy loop indices
624      INTEGER  :: ijpij, i_fill, jl0 
625      REAL(wp) :: zarg, zV, zconv, zdh
626      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
627      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
628      INTEGER , POINTER, DIMENSION(:)         ::   itest
629 
630      CALL wrk_alloc( 4, itest )
631      !--------------------------------------------------------------------
632      ! initialisation of variables
633      !--------------------------------------------------------------------
634      ijpij = SIZE(zhti,1)
635      zht_i(1:ijpij,1:jpl) = 0._wp
636      zht_s(1:ijpij,1:jpl) = 0._wp
637      za_i (1:ijpij,1:jpl) = 0._wp
638
639      ! ----------------------------------------
640      ! distribution over the jpl ice categories
641      ! ----------------------------------------
642      DO ji = 1, ijpij
643         
644         IF( zhti(ji) > 0._wp ) THEN
645
646         ! initialisation of tests
647         itest(:)  = 0
648         
649         i_fill = jpl + 1                                             !====================================
650         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
651            ! iteration                                               !====================================
652            i_fill = i_fill - 1
653           
654            ! initialisation of ice variables for each try
655            zht_i(ji,1:jpl) = 0._wp
656            za_i (ji,1:jpl) = 0._wp
657           
658            ! *** case very thin ice: fill only category 1
659            IF ( i_fill == 1 ) THEN
660               zht_i(ji,1) = zhti(ji)
661               za_i (ji,1) = zai (ji)
662
663            ! *** case ice is thicker: fill categories >1
664            ELSE
665
666               ! Fill ice thicknesses except the last one (i_fill) by hmean
667               DO jl = 1, i_fill - 1
668                  zht_i(ji,jl) = hi_mean(jl)
669               END DO
670               
671               ! find which category (jl0) the input ice thickness falls into
672               jl0 = i_fill
673               DO jl = 1, i_fill
674                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
675                     jl0 = jl
676           CYCLE
677                  ENDIF
678               END DO
679               
680               ! Concentrations in the (i_fill-1) categories
681               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
682               DO jl = 1, i_fill - 1
683                  IF ( jl == jl0 ) CYCLE
684                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
685                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
686               END DO
687               
688               ! Concentration in the last (i_fill) category
689               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
690               
691               ! Ice thickness in the last (i_fill) category
692               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
693               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
694               
695            ENDIF ! case ice is thick or thin
696           
697            !---------------------
698            ! Compatibility tests
699            !---------------------
700            ! Test 1: area conservation
701            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
702            IF ( zconv < epsi06 ) itest(1) = 1
703           
704            ! Test 2: volume conservation
705            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
706            IF ( zconv < epsi06 ) itest(2) = 1
707           
708            ! Test 3: thickness of the last category is in-bounds ?
709            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
710           
711            ! Test 4: positivity of ice concentrations
712            itest(4) = 1
713            DO jl = 1, i_fill
714               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
715            END DO           
716                                                           !============================
717         END DO                                            ! end iteration on categories
718                                                           !============================
719         ENDIF ! if zhti > 0
720      END DO ! i loop
721
722      ! ------------------------------------------------
723      ! Adding Snow in each category where za_i is not 0
724      ! ------------------------------------------------
725      DO jl = 1, jpl
726         DO ji = 1, ijpij
727            IF( za_i(ji,jl) > 0._wp ) THEN
728               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
729               ! In case snow load is in excess that would lead to transformation from snow to ice
730               ! Then, transfer the snow excess into the ice (different from limthd_dh)
731               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
732               ! recompute ht_i, ht_s avoiding out of bounds values
733               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
734               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
735            ENDIF
736         ENDDO
737      ENDDO
738
739      CALL wrk_dealloc( 4, itest )
740      !
741    END SUBROUTINE lim_var_itd
742
743
744#else
745   !!----------------------------------------------------------------------
746   !!   Default option         Dummy module          NO  LIM3 sea-ice model
747   !!----------------------------------------------------------------------
748CONTAINS
749   SUBROUTINE lim_var_agg          ! Empty routines
750   END SUBROUTINE lim_var_agg
751   SUBROUTINE lim_var_glo2eqv      ! Empty routines
752   END SUBROUTINE lim_var_glo2eqv
753   SUBROUTINE lim_var_eqv2glo      ! Empty routines
754   END SUBROUTINE lim_var_eqv2glo
755   SUBROUTINE lim_var_salprof      ! Empty routines
756   END SUBROUTINE lim_var_salprof
757   SUBROUTINE lim_var_bv           ! Emtpy routines
758   END SUBROUTINE lim_var_bv
759   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
760   END SUBROUTINE lim_var_salprof1d
761   SUBROUTINE lim_var_zapsmall
762   END SUBROUTINE lim_var_zapsmall
763   SUBROUTINE lim_var_itd
764   END SUBROUTINE lim_var_itd
765#endif
766
767   !!======================================================================
768END MODULE limvar
Note: See TracBrowser for help on using the repository browser.