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 trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/limvar.F90 @ 1156

Last change on this file since 1156 was 1156, checked in by rblod, 16 years ago

Update Id and licence information, see ticket #210
  • Property svn:keywords set to Id
File size: 26.6 KB
Line 
1MODULE limvar
2   !!----------------------------------------------------------------------
3   !!   'key_lim3'                                      LIM3 sea-ice model
4   !!----------------------------------------------------------------------
5   !!======================================================================
6   !!                       ***  MODULE limvar ***
7   !!                 Different sets of ice model variables
8   !!                   how to switch from one to another
9   !!
10   !!                 There are three sets of variables
11   !!                 VGLO : global variables of the model
12   !!                        - v_i (jpi,jpj,jpl)
13   !!                        - v_s (jpi,jpj,jpl)
14   !!                        - a_i (jpi,jpj,jpl)
15   !!                        - t_s (jpi,jpj,jpl)
16   !!                        - e_i (jpi,jpj,nlay_i,jpl)
17   !!                        - smv_i(jpi,jpj,jpl)
18   !!                        - oa_i (jpi,jpj,jpl)
19   !!                 VEQV : equivalent variables sometimes used in the model
20   !!                        - ht_i(jpi,jpj,jpl)
21   !!                        - ht_s(jpi,jpj,jpl)
22   !!                        - t_i (jpi,jpj,nlay_i,jpl)
23   !!                        ...
24   !!                 VAGG : aggregate variables, averaged/summed over all
25   !!                        thickness categories
26   !!                        - vt_i(jpi,jpj)
27   !!                        - vt_s(jpi,jpj)
28   !!                        - at_i(jpi,jpj)
29   !!                        - et_s(jpi,jpj)  !total snow heat content
30   !!                        - et_i(jpi,jpj)  !total ice thermal content
31   !!                        - smt_i(jpi,jpj) !mean ice salinity
32   !!                        - ot_i(jpi,jpj)  !average ice age
33   !!======================================================================
34#if defined key_lim3
35   !!----------------------------------------------------------------------
36   !! * Modules used
37   USE dom_ice
38   USE par_oce          ! ocean parameters
39   USE phycst           ! physical constants (ocean directory)
40   USE ice_oce          ! ice variables
41   USE sbc_oce          ! Surface boundary condition: ocean fields
42   USE thd_ice
43   USE in_out_manager
44   USE ice
45   USE par_ice
46
47   IMPLICIT NONE
48   PRIVATE
49
50   !! * Routine accessibility
51   PUBLIC lim_var_agg
52   PUBLIC lim_var_glo2eqv
53   PUBLIC lim_var_eqv2glo
54   PUBLIC lim_var_salprof
55   PUBLIC lim_var_bv
56   PUBLIC lim_var_salprof1d
57
58   !! * Module variables
59   REAL(wp)  ::           &  ! constant values
60      epsi20 = 1e-20   ,  &
61      epsi13 = 1e-13   ,  &
62      zzero  = 0.e0    ,  &
63      zone   = 1.e0
64
65   !!----------------------------------------------------------------------
66   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
67   !! $Id$
68   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
69   !!----------------------------------------------------------------------
70
71
72CONTAINS
73
74   SUBROUTINE lim_var_agg(n)
75      !!------------------------------------------------------------------
76      !!                ***  ROUTINE lim_var_agg  ***
77      !! ** Purpose :
78      !!        This routine aggregates ice-thickness-category variables to 
79      !!                                all-ice variables
80      !!        i.e. it turns VGLO into VAGG
81      !! ** Method  :
82      !!
83      !! ** Arguments :
84      !!           kideb , kiut : Starting and ending points on which the
85      !!                         the computation is applied
86      !!
87      !! ** Inputs / Ouputs : (global commons)
88      !! ** Arguments : n = 1, at_i vt_i only
89      !!                n = 2 everything
90      !!
91      !! ** External :
92      !!
93      !! ** References :
94      !!
95      !! ** History :
96      !!           (01-2006) Martin Vancoppenolle, UCL-ASTR
97      !!
98      !! note : you could add an argument when you need only at_i, vt_i
99      !!        and when you need everything
100      !!------------------------------------------------------------------
101      !! * Arguments
102
103      !! * Local variables
104      INTEGER ::   ji,       &   ! spatial dummy loop index
105         jj,       &   ! spatial dummy loop index
106         jk,       &   ! vertical layering dummy loop index
107         jl            ! ice category dummy loop index
108
109      REAL ::      zeps, epsi16, zinda, epsi06
110
111      INTEGER, INTENT( in ) ::   n     ! describes what is needed
112
113      !!-- End of declarations
114      !!----------------------------------------------------------------------------------------------
115      zeps = 1.0e-13
116      epsi16 = 1.0e-16
117      epsi06 = 1.0e-6
118
119      !------------------
120      ! Zero everything
121      !------------------
122
123      vt_i(:,:)  = 0.0
124      vt_s(:,:)  = 0.0
125      at_i(:,:)  = 0.0 
126      ato_i(:,:) = 1.0 
127
128      IF ( n .GT. 1 ) THEN
129         et_s(:,:)  = 0.0
130         ot_i(:,:)  = 0.0
131         smt_i(:,:) = 0.0
132         et_i(:,:)  = 0.0
133      ENDIF
134
135      !--------------------
136      ! Compute variables
137      !--------------------
138
139      DO jl = 1, jpl
140         DO jj = 1, jpj
141            DO ji = 1, jpi
142
143               vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
144               vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
145               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
146
147               zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 
148               icethi(ji,jj) = vt_i(ji,jj) / MAX(at_i(ji,jj),epsi16)*zinda 
149               ! ice thickness
150            END DO
151         END DO
152      END DO
153
154      DO jj = 1, jpj
155         DO ji = 1, jpi
156            ato_i(ji,jj) = MAX(1.0 - at_i(ji,jj), 0.0)   ! open water fraction
157         END DO
158      END DO
159
160      IF ( n .GT. 1 ) THEN
161
162         DO jl = 1, jpl
163            DO jj = 1, jpj
164               DO ji = 1, jpi
165                  et_s(ji,jj)  = et_s(ji,jj)  +     &       ! snow heat content
166                     e_s(ji,jj,1,jl)
167                  zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - 0.10 ) ) 
168                  smt_i(ji,jj) = smt_i(ji,jj) +     &       ! ice salinity
169                     smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , zeps ) * &
170                     zinda
171                  zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) ) 
172                  ot_i(ji,jj)  = ot_i(ji,jj)  +     &       ! ice age
173                     oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , zeps ) * &
174                     zinda
175               END DO
176            END DO
177         END DO
178
179         DO jl = 1, jpl
180            DO jk = 1, nlay_i
181               DO jj = 1, jpj
182                  DO ji = 1, jpi
183                     et_i(ji,jj) = et_i(ji,jj) + e_i(ji,jj,jk,jl) ! ice heat
184                     ! content
185                  END DO
186               END DO
187            END DO
188         END DO
189
190      ENDIF ! n .GT. 1
191
192   END SUBROUTINE lim_var_agg
193
194   !==============================================================================
195
196   SUBROUTINE lim_var_glo2eqv
197      !!------------------------------------------------------------------
198      !!                ***  ROUTINE lim_var_glo2eqv ***'
199      !! ** Purpose :
200      !!        This routine computes equivalent variables as function of   
201      !!                              global variables
202      !!        i.e. it turns VGLO into VEQV
203      !! ** Method  :
204      !!
205      !! ** Arguments :
206      !!           kideb , kiut : Starting and ending points on which the
207      !!                         the computation is applied
208      !!
209      !! ** Inputs / Ouputs :
210      !!
211      !! ** External :
212      !!
213      !! ** References :
214      !!
215      !! ** History :
216      !!           (01-2006) Martin Vancoppenolle, UCL-ASTR
217      !!
218      !!------------------------------------------------------------------
219
220      !! * Local variables
221      INTEGER ::   ji,       &   ! spatial dummy loop index
222         jj,       &   ! spatial dummy loop index
223         jk,       &   ! vertical layering dummy loop index
224         jl            ! ice category dummy loop index
225
226      REAL :: zq_i, zaaa, zbbb, zccc, zdiscrim, &
227         ztmelts, zindb, zq_s, zfac1, zfac2
228
229      REAL :: zeps, epsi06
230
231      zeps    = 1.0e-10
232      epsi06  = 1.0e-06
233
234      !!-- End of declarations
235      !!------------------------------------------------------------------------------
236
237      !-------------------------------------------------------
238      ! Ice thickness, snow thickness, ice salinity, ice age
239      !-------------------------------------------------------
240!CDIR NOVERRCHK
241      DO jl = 1, jpl
242!CDIR NOVERRCHK
243         DO jj = 1, jpj
244!CDIR NOVERRCHK
245            DO ji = 1, jpi
246               zindb          = 1.0-MAX(0.0,SIGN(1.0,- a_i(ji,jj,jl))) !0 if no ice and 1 if yes
247               ht_i(ji,jj,jl) = v_i(ji,jj,jl)   / MAX( a_i(ji,jj,jl) , zeps ) * zindb
248               ht_s(ji,jj,jl) = v_s(ji,jj,jl)   / MAX( a_i(ji,jj,jl) , zeps ) * zindb
249               o_i(ji,jj,jl)  = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , zeps ) * zindb
250            END DO
251         END DO
252      END DO
253
254      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) )THEN
255
256!CDIR NOVERRCHK
257         DO jl = 1, jpl
258!CDIR NOVERRCHK
259            DO jj = 1, jpj
260!CDIR NOVERRCHK
261               DO ji = 1, jpi
262                  zindb          = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes
263                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX(v_i(ji,jj,jl),zeps) * zindb
264               END DO
265            END DO
266         END DO
267
268      ENDIF
269
270      ! salinity profile
271      CALL lim_var_salprof
272
273      !-------------------
274      ! Ice temperatures
275      !-------------------
276!CDIR NOVERRCHK
277      DO jl = 1, jpl
278!CDIR NOVERRCHK
279         DO jk = 1, nlay_i
280!CDIR NOVERRCHK
281            DO jj = 1, jpj
282!CDIR NOVERRCHK
283               DO ji = 1, jpi
284                  !Energy of melting q(S,T) [J.m-3]
285                  zq_i       = e_i(ji,jj,jk,jl) / area(ji,jj) / &
286                     MAX( v_i(ji,jj,jl) , epsi06 ) * nlay_i 
287                  ! zindb = 0 if no ice and 1 if yes
288                  zindb      = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) )
289                  !convert units ! very important that this line is here
290                  zq_i       = zq_i * unit_fac * zindb
291                  !Ice layer melt temperature
292                  ztmelts    =  -tmut*s_i(ji,jj,jk,jl) + rtt
293                  !Conversion q(S,T) -> T (second order equation)
294                  zaaa       =  cpic
295                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + &
296                     zq_i / rhoic - lfus
297                  zccc       =  lfus * (ztmelts-rtt)
298                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) )
299                  t_i(ji,jj,jk,jl) = rtt + zindb *( - zbbb - zdiscrim ) / & 
300                     ( 2.0 *zaaa )
301                  t_i(ji,jj,jk,jl) = MIN( rtt, MAX(173.15, t_i(ji,jj,jk,jl) ) )
302               END DO
303            END DO
304         END DO
305      END DO
306
307      !--------------------
308      ! Snow temperatures
309      !--------------------
310      zfac1 = 1. / ( rhosn * cpic )
311      zfac2 = lfus / cpic 
312!CDIR NOVERRCHK
313      DO jl = 1, jpl
314!CDIR NOVERRCHK
315         DO jk = 1, nlay_s
316!CDIR NOVERRCHK
317            DO jj = 1, jpj
318!CDIR NOVERRCHK
319               DO ji = 1, jpi
320                  !Energy of melting q(S,T) [J.m-3]
321                  zq_s       = e_s(ji,jj,jk,jl) / area(ji,jj) / &
322                     MAX( v_s(ji,jj,jl) , epsi06 ) * nlay_s 
323                  ! zindb = 0 if no ice and 1 if yes
324                  zindb      = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) ) )
325                  !convert units ! very important that this line is here
326                  zq_s       = zq_s * unit_fac * zindb
327                  t_s(ji,jj,jk,jl) = rtt + zindb * ( - zfac1 * zq_s + zfac2 )
328                  t_s(ji,jj,jk,jl) = MIN( rtt, MAX(173.15, t_s(ji,jj,jk,jl) ) )
329
330               END DO
331            END DO
332         END DO
333      END DO
334
335      !-------------------
336      ! Mean temperature
337      !-------------------
338      tm_i(:,:) = 0.0
339!CDIR NOVERRCHK
340      DO jl = 1, jpl
341!CDIR NOVERRCHK
342         DO jk = 1, nlay_i
343!CDIR NOVERRCHK
344            DO jj = 1, jpj
345!CDIR NOVERRCHK
346               DO ji = 1, jpi
347                  zindb          = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)))
348                  zindb          = zindb*1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl)))
349                  tm_i(ji,jj) = tm_i(ji,jj) + t_i(ji,jj,jk,jl)*v_i(ji,jj,jl) / &
350                     REAL(nlay_i) / MAX( vt_i(ji,jj) , zeps )
351               END DO
352            END DO
353         END DO
354      END DO
355
356   END SUBROUTINE lim_var_glo2eqv
357
358   !===============================================================================
359
360   SUBROUTINE lim_var_eqv2glo
361      !!------------------------------------------------------------------
362      !!                ***  ROUTINE lim_var_eqv2glo ***'
363      !! ** Purpose :
364      !!        This routine computes global     variables as function of   
365      !!                              equivalent variables
366      !!        i.e. it turns VEQV into VGLO
367      !! ** Method  :
368      !!
369      !! ** Arguments :
370      !!
371      !! ** Inputs / Ouputs : (global commons)
372      !!
373      !! ** External :
374      !!
375      !! ** References :
376      !!
377      !! ** History :
378      !!           (01-2006) Martin Vancoppenolle, UCL-ASTR
379      !!                     Take it easy man
380      !!                     Life is just a simple game, between
381      !!                     ups / and downs \ :@)
382      !!
383      !!------------------------------------------------------------------
384
385      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
386      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
387      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
388      oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:)
389
390   END SUBROUTINE lim_var_eqv2glo
391
392   !===============================================================================
393
394   SUBROUTINE lim_var_salprof
395      !!------------------------------------------------------------------
396      !!                ***  ROUTINE lim_var_salprof ***'
397      !! ** Purpose :
398      !!        This routine computes salinity profile in function of
399      !!        bulk salinity     
400      !!
401      !! ** Method  : If bulk salinity greater than s_i_1,
402      !!              the profile is assumed to be constant (S_inf)
403      !!              If bulk salinity lower than s_i_0,
404      !!              the profile is linear with 0 at the surface (S_zero)
405      !!              If it is between s_i_0 and s_i_1, it is a
406      !!              alpha-weighted linear combination of s_inf and s_zero
407      !!
408      !! ** References : Vancoppenolle et al., 2007 (in preparation)
409      !!
410      !! ** History :
411      !!           (08-2006) Martin Vancoppenolle, UCL-ASTR
412      !!                     Take it easy man
413      !!                     Life is just a simple game, between ups
414      !!                     / and downs \ :@)
415      !!
416      !!------------------------------------------------------------------
417      !! * Arguments
418
419      !! * Local variables
420      INTEGER ::             &
421         ji            ,     & !: spatial dummy loop index
422         jj            ,     & !: spatial dummy loop index
423         jk            ,     & !: vertical layering dummy loop index
424         jl                    !: ice category dummy loop index
425
426      REAL(wp) ::            &
427         dummy_fac0    ,     & !: dummy factor used in computations
428         dummy_fac1    ,     & !: dummy factor used in computations
429         dummy_fac     ,     & !: dummy factor used in computations
430         zind0         ,     & !: switch, = 1 if sm_i lt s_i_0
431         zind01        ,     & !: switch, = 1 if sm_i between s_i_0 and s_i_1
432         zindbal       ,     & !: switch, = 1, if 2*sm_i gt sss_m
433         zargtemp              !: dummy factor
434
435      REAL(wp), DIMENSION(nlay_i) ::      &
436         zs_zero               !: linear salinity profile for salinities under
437      !: s_i_0
438
439      REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
440         z_slope_s     ,     & !: slope of the salinity profile
441         zalpha                !: weight factor for s between s_i_0 and s_i_1
442
443      !!-- End of declarations
444      !!------------------------------------------------------------------------------
445
446      !---------------------------------------
447      ! Vertically constant, constant in time
448      !---------------------------------------
449
450      IF ( num_sal .EQ. 1 ) THEN
451
452         s_i(:,:,:,:) = bulk_sal
453
454      ENDIF
455
456      !-----------------------------------
457      ! Salinity profile, varying in time
458      !-----------------------------------
459
460      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) )THEN
461
462         DO jk = 1, nlay_i
463            s_i(:,:,jk,:)  = sm_i(:,:,:)
464         END DO ! jk
465
466         ! Slope of the linear profile zs_zero
467         !-------------------------------------
468         DO jl = 1, jpl
469            DO jj = 1, jpj
470               DO ji = 1, jpi
471                  z_slope_s(ji,jj,jl) = 2.0 * sm_i(ji,jj,jl) / MAX( 0.01      &
472                     , ht_i(ji,jj,jl) )
473               END DO ! ji
474            END DO ! jj
475         END DO ! jl
476
477         ! Weighting factor between zs_zero and zs_inf
478         !---------------------------------------------
479         dummy_fac0 = 1. / ( ( s_i_0 - s_i_1 ) )
480         dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )
481
482         zalpha(:,:,:) = 0.0
483
484!CDIR NOVERRCHK
485         DO jl = 1, jpl
486!CDIR NOVERRCHK
487            DO jj = 1, jpj
488!CDIR NOVERRCHK
489               DO ji = 1, jpi
490                  ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
491                  zind0  = MAX( 0.0   , SIGN( 1.0  , s_i_0 - sm_i(ji,jj,jl) ) ) 
492                  ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws
493                  zind01 = ( 1.0 - zind0 ) *                                  &
494                     MAX( 0.0   , SIGN( 1.0  , s_i_1 - sm_i(ji,jj,jl) ) ) 
495                  ! If 2.sm_i GE sss_m then zindbal = 1
496                  zindbal = MAX( 0.0 , SIGN( 1.0 , 2. * sm_i(ji,jj,jl) -      &
497                     sss_m(ji,jj) ) )
498                  zalpha(ji,jj,jl) = zind0  * 1.0                             &
499                     + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + &
500                     dummy_fac1 )
501                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1.0 - zindbal )
502               END DO
503            END DO
504         END DO
505
506         ! Computation of the profile
507         !----------------------------
508         dummy_fac = 1. / nlay_i
509
510         DO jl = 1, jpl
511            DO jk = 1, nlay_i
512               DO jj = 1, jpj
513                  DO ji = 1, jpi
514                     ! linear profile with 0 at the surface
515                     zs_zero(jk)      = z_slope_s(ji,jj,jl) * ( jk - 1./2. ) * &
516                        ht_i(ji,jj,jl) * dummy_fac
517                     ! weighting the profile
518                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero(jk) +       &
519                        ( 1.0 - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
520                  END DO ! ji
521               END DO ! jj
522            END DO ! jk
523         END DO ! jl
524
525      ENDIF ! num_sal
526
527      !-------------------------------------------------------
528      ! Vertically varying salinity profile, constant in time
529      !-------------------------------------------------------
530      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
531
532      IF ( num_sal .EQ. 3 ) THEN
533
534         sm_i(:,:,:) = 2.30
535
536!CDIR NOVERRCHK
537         DO jl = 1, jpl
538!CDIR NOVERRCHK
539            DO jk = 1, nlay_i
540!CDIR NOVERRCHK
541               DO jj = 1, jpj
542!CDIR NOVERRCHK
543                  DO ji = 1, jpi
544                     zargtemp  = ( jk - 0.5 ) / nlay_i
545                     s_i(ji,jj,jk,jl) =  1.6 - 1.6 * COS( 3.14169265 * &
546                        ( zargtemp**(0.407/           &
547                        ( 0.573 + zargtemp ) ) ) )
548                  END DO ! ji
549               END DO ! jj
550            END DO ! jk
551         END DO ! jl
552
553      ENDIF ! num_sal
554
555   END SUBROUTINE lim_var_salprof
556
557   !===============================================================================
558
559   SUBROUTINE lim_var_bv
560      !!------------------------------------------------------------------
561      !!                ***  ROUTINE lim_var_bv ***'
562      !! ** Purpose :
563      !!        This routine computes mean brine volume (%) in sea ice
564      !!
565      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
566      !!
567      !! ** Arguments :
568      !!
569      !! ** Inputs / Ouputs : (global commons)
570      !!
571      !! ** External :
572      !!
573      !! ** References : Vancoppenolle et al., JGR, 2007
574      !!
575      !! ** History :
576      !!           (08-2006) Martin Vancoppenolle, UCL-ASTR
577      !!
578      !!------------------------------------------------------------------
579      !! * Arguments
580
581      !! * Local variables
582      INTEGER ::   ji,       &   ! spatial dummy loop index
583         jj,       &   ! spatial dummy loop index
584         jk,       &   ! vertical layering dummy loop index
585         jl            ! ice category dummy loop index
586
587      REAL :: zbvi,          &   ! brine volume for a single ice category
588         zeps,          &   ! very small value
589         zindb              ! is there ice or not
590
591      !!-- End of declarations
592      !!------------------------------------------------------------------------------
593
594      zeps = 1.0e-13
595      bv_i(:,:) = 0.0
596!CDIR NOVERRCHK
597      DO jl = 1, jpl
598!CDIR NOVERRCHK
599         DO jk = 1, nlay_i
600!CDIR NOVERRCHK
601            DO jj = 1, jpj
602!CDIR NOVERRCHK
603               DO ji = 1, jpi
604                  zindb          = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes
605                  zbvi = - zindb * tmut *s_i(ji,jj,jk,jl) /             & 
606                     MIN( t_i(ji,jj,jk,jl) - 273.15 , zeps )         &
607                     * v_i(ji,jj,jl) / REAL(nlay_i)
608                  bv_i(ji,jj) = bv_i(ji,jj) + zbvi  &
609                     / MAX( vt_i(ji,jj) , zeps )
610               END DO
611            END DO
612         END DO
613      END DO
614
615   END SUBROUTINE lim_var_bv
616
617   !===============================================================================
618
619   SUBROUTINE lim_var_salprof1d(kideb,kiut)
620      !!-------------------------------------------------------------------
621      !!                  ***  ROUTINE lim_thd_salprof1d  ***
622      !!
623      !! ** Purpose :   1d computation of the sea ice salinity profile
624      !!                Works with 1d vectors and is used by thermodynamic
625      !!                modules
626      !!
627      !! history :
628      !!   3.0  !  May  2007 M. Vancoppenolle  Original code
629      !!-------------------------------------------------------------------
630      INTEGER, INTENT(in) :: &
631         kideb, kiut             ! thickness category index
632
633      INTEGER ::             &
634         ji, jk,             &   ! geographic and layer index
635         zji, zjj
636
637      REAL(wp) ::            &
638         dummy_fac0,         &   ! dummy factors
639         dummy_fac1,         &
640         dummy_fac2,         &
641         zalpha    ,         &   ! weighting factor
642         zind0     ,         &   ! switches as in limvar
643         zind01    ,         &   ! switch
644         zindbal   ,         &   ! switch if in freshwater area
645         zargtemp
646
647      REAL(wp), DIMENSION(jpij) ::            &
648         z_slope_s
649
650      REAL(wp), DIMENSION(jpij,jkmax) ::      &
651         zs_zero
652      !!-------------------------------------------------------------------
653
654      !---------------------------------------
655      ! Vertically constant, constant in time
656      !---------------------------------------
657
658      IF ( num_sal .EQ. 1 ) THEN
659
660         s_i_b(:,:) = bulk_sal
661
662      ENDIF
663
664      !------------------------------------------------------
665      ! Vertically varying salinity profile, varying in time
666      !------------------------------------------------------
667
668      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
669
670         ! Slope of the linear profile zs_zero
671         !-------------------------------------
672!CDIR NOVERRCHK
673         DO ji = kideb, kiut 
674            z_slope_s(ji) = 2.0 * sm_i_b(ji) / MAX( 0.01      &
675               , ht_i_b(ji) )
676         END DO ! ji
677
678         ! Weighting factor between zs_zero and zs_inf
679         !---------------------------------------------
680         dummy_fac0 = 1. / ( ( s_i_0 - s_i_1 ) )
681         dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )
682         dummy_fac2 = 1. / nlay_i
683
684!CDIR NOVERRCHK
685         DO jk = 1, nlay_i
686!CDIR NOVERRCHK
687            DO ji = kideb, kiut
688               zji    =  MOD( npb(ji) - 1, jpi ) + 1
689               zjj    =  ( npb(ji) - 1 ) / jpi + 1
690               zalpha = 0.0
691               ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
692               zind0  = MAX( 0.0   , SIGN( 1.0  , s_i_0 - sm_i_b(ji) ) ) 
693               ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws
694               zind01 = ( 1.0 - zind0 ) *                                  &
695                  MAX( 0.0   , SIGN( 1.0  , s_i_1 - sm_i_b(ji) ) ) 
696               ! if 2.sm_i GE sss_m then zindbal = 1
697               zindbal = MAX( 0.0 , SIGN( 1.0 , 2. * sm_i_b(ji) -      &
698                  sss_m(zji,zjj) ) )
699
700               zalpha = zind0  * 1.0                                       &
701                  + zind01 * ( sm_i_b(ji) * dummy_fac0 +           &
702                  dummy_fac1 )
703               zalpha = zalpha * ( 1.0 - zindbal )
704
705               zs_zero(ji,jk) = z_slope_s(ji) * ( jk - 1./2. ) * &
706                  ht_i_b(ji) * dummy_fac2
707               ! weighting the profile
708               s_i_b(ji,jk) = zalpha * zs_zero(ji,jk) +       &
709                  ( 1.0 - zalpha ) * sm_i_b(ji)
710            END DO ! ji
711         END DO ! jk
712
713      ENDIF ! num_sal
714
715      !-------------------------------------------------------
716      ! Vertically varying salinity profile, constant in time
717      !-------------------------------------------------------
718      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
719
720      IF ( num_sal .EQ. 3 ) THEN
721
722         sm_i_b(:) = 2.30
723
724!CDIR NOVERRCHK
725         DO ji = kideb, kiut
726!CDIR NOVERRCHK
727            DO jk = 1, nlay_i
728               zargtemp  = ( jk - 0.5 ) / nlay_i
729               s_i_b(ji,jk)  =  1.6 - 1.6*cos(3.14169265*(zargtemp**(0.407/ &
730                  (0.573+zargtemp))))
731            END DO ! jk
732         END DO ! ji
733
734      ENDIF ! num_sal
735
736   END SUBROUTINE lim_var_salprof1d
737
738   !===============================================================================
739
740#else
741   !!======================================================================
742   !!                       ***  MODULE limvar  ***
743   !!                          no sea ice model
744   !!======================================================================
745CONTAINS
746   SUBROUTINE lim_var_agg          ! Empty routines
747   END SUBROUTINE lim_var_agg
748   SUBROUTINE lim_var_glo2eqv      ! Empty routines
749   END SUBROUTINE lim_var_glo2eqv
750   SUBROUTINE lim_var_eqv2glo      ! Empty routines
751   END SUBROUTINE lim_var_eqv2glo
752   SUBROUTINE lim_var_salprof      ! Empty routines
753   END SUBROUTINE lim_var_salprof
754   SUBROUTINE lim_var_bv           ! Emtpy routines
755   END SUBROUTINE lim_var_bv
756   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
757   END SUBROUTINE lim_var_salprof1d
758
759#endif
760END MODULE limvar
Note: See TracBrowser for help on using the repository browser.