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

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 7881

Last change on this file since 7881 was 7813, checked in by clem, 7 years ago

synchronize trunk with 3.6

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