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

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

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