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

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

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