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

source: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 7773

Last change on this file since 7773 was 7773, checked in by mattmartin, 7 years ago

Committing updates after doing the following:

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