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

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limvar.F90 @ 830

Last change on this file since 830 was 825, checked in by ctlod, 16 years ago

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