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/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 3963

Last change on this file since 3963 was 3963, checked in by clem, 11 years ago

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116
  • Property svn:keywords set to Id
File size: 23.0 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   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
33   !!----------------------------------------------------------------------
34#if defined key_lim3
35   !!----------------------------------------------------------------------
36   !!   'key_lim3'                                      LIM3 sea-ice model
37   !!----------------------------------------------------------------------
38   !!   lim_var_agg       :
39   !!   lim_var_glo2eqv   :
40   !!   lim_var_eqv2glo   :
41   !!   lim_var_salprof   :
42   !!   lim_var_salprof1d :
43   !!   lim_var_bv        :
44   !!----------------------------------------------------------------------
45   USE par_oce          ! ocean parameters
46   USE phycst           ! physical constants (ocean directory)
47   USE sbc_oce          ! Surface boundary condition: ocean fields
48   USE ice              ! LIM variables
49   USE par_ice          ! LIM parameters
50   USE dom_ice          ! LIM domain
51   USE thd_ice          ! LIM thermodynamics
52   USE in_out_manager   ! I/O manager
53   USE lib_mpp          ! MPP library
54   USE wrk_nemo         ! work arrays
55   USE lib_fortran      ! to use key_nosignedzero
56
57   IMPLICIT NONE
58   PRIVATE
59
60   PUBLIC   lim_var_agg          !
61   PUBLIC   lim_var_glo2eqv      !
62   PUBLIC   lim_var_eqv2glo      !
63   PUBLIC   lim_var_salprof      !
64   PUBLIC   lim_var_icetm        !
65   PUBLIC   lim_var_bv           !
66   PUBLIC   lim_var_salprof1d    !
67
68   REAL(wp) ::   epsi20 = 1.e-20_wp   ! module constants
69   REAL(wp) ::   epsi16 = 1.e-16_wp   !    -       -
70   REAL(wp) ::   epsi13 = 1.e-13_wp   !    -       -
71   REAL(wp) ::   epsi10 = 1.e-10_wp   !    -       -
72   REAL(wp) ::   epsi06 = 1.e-06_wp   !    -       -
73   REAL(wp) ::   zzero = 0.e0        !    -       -
74   REAL(wp) ::   zone  = 1.e0        !    -       -
75
76   !!----------------------------------------------------------------------
77   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
78   !! $Id$
79   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
80   !!----------------------------------------------------------------------
81CONTAINS
82
83   SUBROUTINE lim_var_agg( kn )
84      !!------------------------------------------------------------------
85      !!                ***  ROUTINE lim_var_agg  ***
86      !!
87      !! ** Purpose :   aggregates ice-thickness-category variables to all-ice variables
88      !!              i.e. it turns VGLO into VAGG
89      !! ** Method  :
90      !!
91      !! ** Arguments : n = 1, at_i vt_i only
92      !!                n = 2 everything
93      !!
94      !! note : you could add an argument when you need only at_i, vt_i
95      !!        and when you need everything
96      !!------------------------------------------------------------------
97      INTEGER, INTENT( in ) ::   kn     ! =1 at_i & vt only ; = what is needed
98      !
99      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
100      REAL(wp) ::   zinda, zindb
101      !!------------------------------------------------------------------
102
103      !--------------------
104      ! Compute variables
105      !--------------------
106      vt_i (:,:) = 0._wp
107      vt_s (:,:) = 0._wp
108      at_i (:,:) = 0._wp
109      ato_i(:,:) = 1._wp
110      !
111      DO jl = 1, jpl
112         DO jj = 1, jpj
113            DO ji = 1, jpi
114               !
115               vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
116               vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
117               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
118               !
119               zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) ) 
120               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness
121            END DO
122         END DO
123      END DO
124
125      DO jj = 1, jpj
126         DO ji = 1, jpi
127            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   ! open water fraction
128         END DO
129      END DO
130
131      IF( kn > 1 ) THEN
132         et_s (:,:) = 0._wp
133         ot_i (:,:) = 0._wp
134         smt_i(:,:) = 0._wp
135         et_i (:,:) = 0._wp
136         !
137         DO jl = 1, jpl
138            DO jj = 1, jpj
139               DO ji = 1, jpi
140                  zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi16 ) ) 
141                  zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) ) 
142                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content
143                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi16 ) * zinda   ! ice salinity
144                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi16 ) * zindb   ! ice age
145               END DO
146            END DO
147         END DO
148          !
149         DO jl = 1, jpl
150            DO jk = 1, nlay_i
151               et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl)       ! ice heat content
152            END DO
153         END DO
154         !
155      ENDIF
156      !
157   END SUBROUTINE lim_var_agg
158
159
160   SUBROUTINE lim_var_glo2eqv
161      !!------------------------------------------------------------------
162      !!                ***  ROUTINE lim_var_glo2eqv ***
163      !!
164      !! ** Purpose :   computes equivalent variables as function of global variables
165      !!              i.e. it turns VGLO into VEQV
166      !!------------------------------------------------------------------
167      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
168      REAL(wp) ::   zq_i, zaaa, zbbb, zccc, zdiscrim     ! local scalars
169      REAL(wp) ::   ztmelts, zindb, zq_s, zfac1, zfac2   !   -      -
170      !!------------------------------------------------------------------
171
172      !-------------------------------------------------------
173      ! Ice thickness, snow thickness, ice salinity, ice age
174      !-------------------------------------------------------
175      DO jl = 1, jpl
176         DO jj = 1, jpj
177            DO ji = 1, jpi
178               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes
179               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
180               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
181               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
182               a_i(ji,jj,jl) = a_i (ji,jj,jl) * zindb ! clem correction
183            END DO
184         END DO
185      END DO
186
187      IF(  num_sal == 2  .OR.  num_sal == 4  )THEN
188         DO jl = 1, jpl
189            DO jj = 1, jpj
190               DO ji = 1, jpi
191                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes
192                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * zindb
193               END DO
194            END DO
195         END DO
196      ENDIF
197
198      CALL lim_var_salprof      ! salinity profile
199
200      !-------------------
201      ! Ice temperatures
202      !-------------------
203!CDIR NOVERRCHK
204      DO jl = 1, jpl
205!CDIR NOVERRCHK
206         DO jk = 1, nlay_i
207!CDIR NOVERRCHK
208            DO jj = 1, jpj
209!CDIR NOVERRCHK
210               DO ji = 1, jpi
211                  !                                                              ! Energy of melting q(S,T) [J.m-3]
212                  zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi06 ) * REAL(nlay_i,wp) 
213                  zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes
214                  zq_i    = zq_i * unit_fac * zindb                              !convert units
215                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature
216                  !
217                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
218                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i / rhoic - lfus
219                  zccc       =  lfus * (ztmelts-rtt)
220                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
221                  t_i(ji,jj,jk,jl) = rtt + zindb *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )
222                  t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rtt < t_i < rtt
223               END DO
224            END DO
225         END DO
226      END DO
227
228      !--------------------
229      ! Snow temperatures
230      !--------------------
231      zfac1 = 1._wp / ( rhosn * cpic )
232      zfac2 = lfus / cpic 
233      DO jl = 1, jpl
234         DO jk = 1, nlay_s
235            DO jj = 1, jpj
236               DO ji = 1, jpi
237                  !Energy of melting q(S,T) [J.m-3]
238                  zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL(nlay_s,wp)
239                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes
240                  zq_s  = zq_s * unit_fac * zindb                                    ! convert units
241                  !
242                  t_s(ji,jj,jk,jl) = rtt + zindb * ( - zfac1 * zq_s + zfac2 )
243                  t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rtt < t_i < rtt
244               END DO
245            END DO
246         END DO
247      END DO
248
249      !-------------------
250      ! Mean temperature
251      !-------------------
252      tm_i(:,:) = 0._wp
253      DO jl = 1, jpl
254         DO jk = 1, nlay_i
255            DO jj = 1, jpj
256               DO ji = 1, jpi
257                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
258                  tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
259                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
260               END DO
261            END DO
262         END DO
263      END DO
264      !
265   END SUBROUTINE lim_var_glo2eqv
266
267
268   SUBROUTINE lim_var_eqv2glo
269      !!------------------------------------------------------------------
270      !!                ***  ROUTINE lim_var_eqv2glo ***
271      !!
272      !! ** Purpose :   computes global variables as function of equivalent variables
273      !!                i.e. it turns VEQV into VGLO
274      !! ** Method  :
275      !!
276      !! ** History :  (01-2006) Martin Vancoppenolle, UCL-ASTR
277      !!------------------------------------------------------------------
278      !
279      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
280      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
281      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
282      oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:)
283      !
284   END SUBROUTINE lim_var_eqv2glo
285
286
287   SUBROUTINE lim_var_salprof
288      !!------------------------------------------------------------------
289      !!                ***  ROUTINE lim_var_salprof ***
290      !!
291      !! ** Purpose :   computes salinity profile in function of bulk salinity     
292      !!
293      !! ** Method  : If bulk salinity greater than s_i_1,
294      !!              the profile is assumed to be constant (S_inf)
295      !!              If bulk salinity lower than s_i_0,
296      !!              the profile is linear with 0 at the surface (S_zero)
297      !!              If it is between s_i_0 and s_i_1, it is a
298      !!              alpha-weighted linear combination of s_inf and s_zero
299      !!
300      !! ** References : Vancoppenolle et al., 2007 (in preparation)
301      !!------------------------------------------------------------------
302      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
303      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac, zsal      ! local scalar
304      REAL(wp) ::   zind0, zind01, zindbal, zargtemp , zs_zero   !   -      -
305      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha   ! 3D pointer
306      !!------------------------------------------------------------------
307
308      CALL wrk_alloc( jpi, jpj, jpl, z_slope_s, zalpha )
309
310      !---------------------------------------
311      ! Vertically constant, constant in time
312      !---------------------------------------
313      IF( num_sal == 1 )   s_i(:,:,:,:) = bulk_sal
314
315      !-----------------------------------
316      ! Salinity profile, varying in time
317      !-----------------------------------
318
319      IF(   num_sal == 2  .OR.   num_sal == 4   ) THEN
320         !
321         DO jk = 1, nlay_i
322            s_i(:,:,jk,:)  = sm_i(:,:,:)
323         END DO
324         !
325         DO jl = 1, jpl                               ! Slope of the linear profile
326            DO jj = 1, jpj
327               DO ji = 1, jpi
328                  z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( 0.01 , ht_i(ji,jj,jl) )
329               END DO
330            END DO
331         END DO
332         !
333         dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 )       ! Weighting factor between zs_zero and zs_inf
334         dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )
335
336         zalpha(:,:,:) = 0._wp
337         DO jl = 1, jpl
338            DO jj = 1, jpj
339               DO ji = 1, jpi
340                  ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
341                  zind0  = MAX( 0._wp   , SIGN( 1._wp  , s_i_0 - sm_i(ji,jj,jl) ) ) 
342                  ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws
343                  zind01 = ( 1._wp - zind0 ) * MAX( 0._wp   , SIGN( 1._wp  , s_i_1 - sm_i(ji,jj,jl) ) ) 
344                  ! If 2.sm_i GE sss_m then zindbal = 1
345                  zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
346                  zalpha(ji,jj,jl) = zind0  + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )
347                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zindbal )
348               END DO
349            END DO
350         END DO
351
352         dummy_fac = 1._wp / REAL( nlay_i )                   ! Computation of the profile
353         DO jl = 1, jpl
354            DO jk = 1, nlay_i
355               DO jj = 1, jpj
356                  DO ji = 1, jpi
357                     !                                      ! linear profile with 0 at the surface
358                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac
359                     !                                      ! weighting the profile
360                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
361                  END DO ! ji
362               END DO ! jj
363            END DO ! jk
364         END DO ! jl
365
366      ENDIF ! num_sal
367
368      !-------------------------------------------------------
369      ! Vertically varying salinity profile, constant in time
370      !-------------------------------------------------------
371
372      IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
373         !
374         sm_i(:,:,:) = 2.30_wp
375         !
376         DO jl = 1, jpl
377!CDIR NOVERRCHK
378            DO jk = 1, nlay_i
379               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
380               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
381               s_i(:,:,jk,jl) =  zsal
382            END DO
383         END DO
384
385      ENDIF ! num_sal
386      !
387      CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
388      !
389   END SUBROUTINE lim_var_salprof
390
391
392   SUBROUTINE lim_var_icetm
393      !!------------------------------------------------------------------
394      !!                ***  ROUTINE lim_var_icetm ***
395      !!
396      !! ** Purpose :   computes mean sea ice temperature
397      !!------------------------------------------------------------------
398      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
399      REAL(wp) ::   zindb   !   -      -
400      !!------------------------------------------------------------------
401
402      ! Mean sea ice temperature
403      tm_i(:,:) = 0._wp
404      DO jl = 1, jpl
405         DO jk = 1, nlay_i
406            DO jj = 1, jpj
407               DO ji = 1, jpi
408                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
409                  tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
410                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
411               END DO
412            END DO
413         END DO
414      END DO
415
416   END SUBROUTINE lim_var_icetm
417
418
419   SUBROUTINE lim_var_bv
420      !!------------------------------------------------------------------
421      !!                ***  ROUTINE lim_var_bv ***
422      !!
423      !! ** Purpose :   computes mean brine volume (%) in sea ice
424      !!
425      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
426      !!
427      !! References : Vancoppenolle et al., JGR, 2007
428      !!------------------------------------------------------------------
429      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
430      REAL(wp) ::   zbvi, zinda, zindb      ! local scalars
431      !!------------------------------------------------------------------
432      !
433      bv_i(:,:) = 0._wp
434      DO jl = 1, jpl
435         DO jk = 1, nlay_i
436            DO jj = 1, jpj
437               DO ji = 1, jpi
438                  zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi16 ) )  )
439                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi16 ) )  )
440                  zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi16 )   &
441                     &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp)
442                  bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi16 )
443               END DO
444            END DO
445         END DO
446      END DO
447      !
448   END SUBROUTINE lim_var_bv
449
450
451   SUBROUTINE lim_var_salprof1d( kideb, kiut )
452      !!-------------------------------------------------------------------
453      !!                  ***  ROUTINE lim_thd_salprof1d  ***
454      !!
455      !! ** Purpose :   1d computation of the sea ice salinity profile
456      !!                Works with 1d vectors and is used by thermodynamic modules
457      !!-------------------------------------------------------------------
458      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
459      !
460      INTEGER  ::   ji, jk    ! dummy loop indices
461      INTEGER  ::   zji, zjj  ! local integers
462      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars
463      REAL(wp) ::   zalpha, zind0, zind01, zindbal, zs_zero              !   -      -
464      !
465      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
466      !!---------------------------------------------------------------------
467
468      CALL wrk_alloc( jpij, z_slope_s )
469
470      !---------------------------------------
471      ! Vertically constant, constant in time
472      !---------------------------------------
473      IF( num_sal == 1 )   s_i_b(:,:) = bulk_sal
474
475      !------------------------------------------------------
476      ! Vertically varying salinity profile, varying in time
477      !------------------------------------------------------
478
479      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
480         !
481         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
482            z_slope_s(ji) = 2._wp * sm_i_b(ji) / MAX( 0.01 , ht_i_b(ji) )
483         END DO
484
485         ! Weighting factor between zs_zero and zs_inf
486         !---------------------------------------------
487         dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 )
488         dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )
489         dummy_fac2 = 1._wp / REAL(nlay_i,wp)
490
491!CDIR NOVERRCHK
492         DO jk = 1, nlay_i
493!CDIR NOVERRCHK
494            DO ji = kideb, kiut
495               zji =  MOD( npb(ji) - 1 , jpi ) + 1
496               zjj =     ( npb(ji) - 1 ) / jpi + 1
497               ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
498               zind0  = MAX( 0._wp , SIGN( 1._wp  , s_i_0 - sm_i_b(ji) ) ) 
499               ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws
500               zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) ) 
501               ! if 2.sm_i GE sss_m then zindbal = 1
502               zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(zji,zjj) ) )
503               !
504               zalpha = (  zind0 + zind01 * ( sm_i_b(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zindbal )
505               !
506               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_b(ji) * dummy_fac2
507               ! weighting the profile
508               s_i_b(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_b(ji)
509            END DO ! ji
510         END DO ! jk
511
512      ENDIF ! num_sal
513
514      !-------------------------------------------------------
515      ! Vertically varying salinity profile, constant in time
516      !-------------------------------------------------------
517
518      IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
519         !
520         sm_i_b(:) = 2.30_wp
521         !
522!CDIR NOVERRCHK
523         DO jk = 1, nlay_i
524            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
525            zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
526            DO ji = kideb, kiut
527               s_i_b(ji,jk) = zsal
528            END DO
529         END DO
530         !
531      ENDIF
532      !
533      CALL wrk_dealloc( jpij, z_slope_s )
534      !
535   END SUBROUTINE lim_var_salprof1d
536
537#else
538   !!----------------------------------------------------------------------
539   !!   Default option         Dummy module          NO  LIM3 sea-ice model
540   !!----------------------------------------------------------------------
541CONTAINS
542   SUBROUTINE lim_var_agg          ! Empty routines
543   END SUBROUTINE lim_var_agg
544   SUBROUTINE lim_var_glo2eqv      ! Empty routines
545   END SUBROUTINE lim_var_glo2eqv
546   SUBROUTINE lim_var_eqv2glo      ! Empty routines
547   END SUBROUTINE lim_var_eqv2glo
548   SUBROUTINE lim_var_salprof      ! Empty routines
549   END SUBROUTINE lim_var_salprof
550   SUBROUTINE lim_var_bv           ! Emtpy routines
551   END SUBROUTINE lim_var_bv
552   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
553   END SUBROUTINE lim_var_salprof1d
554#endif
555
556   !!======================================================================
557END MODULE limvar
Note: See TracBrowser for help on using the repository browser.