source: branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 5350

Last change on this file since 5350 was 5350, checked in by hadcv, 5 years ago

Update to head of the trunk (r5344).

  • Property svn:keywords set to Id
File size: 34.1 KB
Line 
1MODULE limvar
2   !!======================================================================
3   !!                       ***  MODULE limvar ***
4   !!                 Different sets of ice model variables
5   !!                   how to switch from one to another
6   !!
7   !!                 There are three sets of variables
8   !!                 VGLO : global variables of the model
9   !!                        - v_i (jpi,jpj,jpl)
10   !!                        - v_s (jpi,jpj,jpl)
11   !!                        - a_i (jpi,jpj,jpl)
12   !!                        - t_s (jpi,jpj,jpl)
13   !!                        - e_i (jpi,jpj,nlay_i,jpl)
14   !!                        - smv_i(jpi,jpj,jpl)
15   !!                        - oa_i (jpi,jpj,jpl)
16   !!                 VEQV : equivalent variables sometimes used in the model
17   !!                        - ht_i(jpi,jpj,jpl)
18   !!                        - ht_s(jpi,jpj,jpl)
19   !!                        - t_i (jpi,jpj,nlay_i,jpl)
20   !!                        ...
21   !!                 VAGG : aggregate variables, averaged/summed over all
22   !!                        thickness categories
23   !!                        - vt_i(jpi,jpj)
24   !!                        - vt_s(jpi,jpj)
25   !!                        - at_i(jpi,jpj)
26   !!                        - et_s(jpi,jpj)  !total snow heat content
27   !!                        - et_i(jpi,jpj)  !total ice thermal content
28   !!                        - smt_i(jpi,jpj) !mean ice salinity
29   !!                        - ot_i(jpi,jpj)  !average ice age
30   !!======================================================================
31   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
32   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
33   !!----------------------------------------------------------------------
34#if defined key_lim3
35   !!----------------------------------------------------------------------
36   !!   'key_lim3'                                      LIM3 sea-ice model
37   !!----------------------------------------------------------------------
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) 
48
49   IMPLICIT NONE
50   PRIVATE
51
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
61
62   !!----------------------------------------------------------------------
63   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
64   !! $Id$
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE lim_var_agg( kn )
70      !!------------------------------------------------------------------
71      !!                ***  ROUTINE lim_var_agg  ***
72      !!
73      !! ** Purpose :   aggregates ice-thickness-category variables to all-ice variables
74      !!              i.e. it turns VGLO into VAGG
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      !!------------------------------------------------------------------
83      INTEGER, INTENT( in ) ::   kn     ! =1 at_i & vt only ; = what is needed
84      !
85      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
86      !!------------------------------------------------------------------
87
88      !--------------------
89      ! Compute variables
90      !--------------------
91      vt_i (:,:) = 0._wp
92      vt_s (:,:) = 0._wp
93      at_i (:,:) = 0._wp
94      ato_i(:,:) = 1._wp
95      !
96      DO jl = 1, jpl
97         DO jj = 1, jpj
98            DO ji = 1, jpi
99               !
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
103               !
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
106            END DO
107         END DO
108      END DO
109
110      DO jj = 1, jpj
111         DO ji = 1, jpi
112            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   ! open water fraction
113         END DO
114      END DO
115
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         !
122         DO jl = 1, jpl
123            DO jj = 1, jpj
124               DO ji = 1, jpi
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
130               END DO
131            END DO
132         END DO
133         !
134         DO jl = 1, jpl
135            DO jk = 1, nlay_i
136               et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl)       ! ice heat content
137            END DO
138         END DO
139         !
140      ENDIF
141      !
142   END SUBROUTINE lim_var_agg
143
144
145   SUBROUTINE lim_var_glo2eqv
146      !!------------------------------------------------------------------
147      !!                ***  ROUTINE lim_var_glo2eqv ***
148      !!
149      !! ** Purpose :   computes equivalent variables as function of global variables
150      !!              i.e. it turns VGLO into VEQV
151      !!------------------------------------------------------------------
152      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
153      REAL(wp) ::   zq_i, zaaa, zbbb, zccc, zdiscrim     ! local scalars
154      REAL(wp) ::   ztmelts, zq_s, zfac1, zfac2   !   -      -
155      !!------------------------------------------------------------------
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
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
167            END DO
168         END DO
169      END DO
170
171      IF(  nn_icesal == 2  )THEN
172         DO jl = 1, jpl
173            DO jj = 1, jpj
174               DO ji = 1, jpi
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 )
179               END DO
180            END DO
181         END DO
182      ENDIF
183
184      CALL lim_var_salprof      ! salinity profile
185
186      !-------------------
187      ! Ice temperatures
188      !-------------------
189      DO jl = 1, jpl
190         DO jk = 1, nlay_i
191            DO jj = 1, jpj
192               DO ji = 1, jpi
193                  !                                                              ! Energy of melting q(S,T) [J.m-3]
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
197                  !
198                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
199                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus
200                  zccc       =  lfus * (ztmelts-rt0)
201                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
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
204               END DO
205            END DO
206         END DO
207      END DO
208
209      !--------------------
210      ! Snow temperatures
211      !--------------------
212      zfac1 = 1._wp / ( rhosn * cpic )
213      zfac2 = lfus / cpic 
214      DO jl = 1, jpl
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]
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)
221                  !
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
224               END DO
225            END DO
226         END DO
227      END DO
228
229      !-------------------
230      ! Mean temperature
231      !-------------------
232      vt_i (:,:) = 0._wp
233      DO jl = 1, jpl
234         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
235      END DO
236
237      tm_i(:,:) = 0._wp
238      DO jl = 1, jpl
239         DO jk = 1, nlay_i
240            DO jj = 1, jpj
241               DO ji = 1, jpi
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 )
245               END DO
246            END DO
247         END DO
248      END DO
249      tm_i = tm_i + rt0
250      !
251   END SUBROUTINE lim_var_glo2eqv
252
253
254   SUBROUTINE lim_var_eqv2glo
255      !!------------------------------------------------------------------
256      !!                ***  ROUTINE lim_var_eqv2glo ***
257      !!
258      !! ** Purpose :   computes global variables as function of equivalent variables
259      !!                i.e. it turns VEQV into VGLO
260      !! ** Method  :
261      !!
262      !! ** History :  (01-2006) Martin Vancoppenolle, UCL-ASTR
263      !!------------------------------------------------------------------
264      !
265      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
266      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
267      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
268      !
269   END SUBROUTINE lim_var_eqv2glo
270
271
272   SUBROUTINE lim_var_salprof
273      !!------------------------------------------------------------------
274      !!                ***  ROUTINE lim_var_salprof ***
275      !!
276      !! ** Purpose :   computes salinity profile in function of bulk salinity     
277      !!
278      !! ** Method  : If bulk salinity greater than zsi1,
279      !!              the profile is assumed to be constant (S_inf)
280      !!              If bulk salinity lower than zsi0,
281      !!              the profile is linear with 0 at the surface (S_zero)
282      !!              If it is between zsi0 and zsi1, it is a
283      !!              alpha-weighted linear combination of s_inf and s_zero
284      !!
285      !! ** References : Vancoppenolle et al., 2007
286      !!------------------------------------------------------------------
287      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
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
293      !!------------------------------------------------------------------
294
295      CALL wrk_alloc( jpi, jpj, jpl, z_slope_s, zalpha )
296
297      !---------------------------------------
298      ! Vertically constant, constant in time
299      !---------------------------------------
300      IF(  nn_icesal == 1  )   s_i(:,:,:,:) = rn_icesal
301
302      !-----------------------------------
303      ! Salinity profile, varying in time
304      !-----------------------------------
305      IF(  nn_icesal == 2  ) THEN
306         !
307         DO jk = 1, nlay_i
308            s_i(:,:,jk,:)  = sm_i(:,:,:)
309         END DO
310         !
311         DO jl = 1, jpl                               ! Slope of the linear profile
312            DO jj = 1, jpj
313               DO ji = 1, jpi
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) )
316               END DO
317            END DO
318         END DO
319         !
320         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
321         zfac1 = zsi1  / ( zsi1 - zsi0 )
322         !
323         zalpha(:,:,:) = 0._wp
324         DO jl = 1, jpl
325            DO jj = 1, jpj
326               DO ji = 1, jpi
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
332                  ! this is to force a constant salinity profile in the Baltic Sea
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 )
336               END DO
337            END DO
338         END DO
339
340         ! Computation of the profile
341         DO jl = 1, jpl
342            DO jk = 1, nlay_i
343               DO jj = 1, jpj
344                  DO ji = 1, jpi
345                     !                                      ! linear profile with 0 at the surface
346                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
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)
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
355         !
356      ENDIF ! nn_icesal
357
358      !-------------------------------------------------------
359      ! Vertically varying salinity profile, constant in time
360      !-------------------------------------------------------
361
362      IF(  nn_icesal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
363         !
364         sm_i(:,:,:) = 2.30_wp
365         !
366         DO jl = 1, jpl
367            DO jk = 1, nlay_i
368               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
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
373         !
374      ENDIF ! nn_icesal
375      !
376      CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
377      !
378   END SUBROUTINE lim_var_salprof
379
380
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
391      vt_i (:,:) = 0._wp
392      DO jl = 1, jpl
393         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
394      END DO
395
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
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 )
404               END DO
405            END DO
406         END DO
407      END DO
408      tm_i = tm_i + rt0
409
410   END SUBROUTINE lim_var_icetm
411
412
413   SUBROUTINE lim_var_bv
414      !!------------------------------------------------------------------
415      !!                ***  ROUTINE lim_var_bv ***
416      !!
417      !! ** Purpose :   computes mean brine volume (%) in sea ice
418      !!
419      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
420      !!
421      !! References : Vancoppenolle et al., JGR, 2007
422      !!------------------------------------------------------------------
423      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
424      REAL(wp) ::   zbvi             ! local scalars
425      !!------------------------------------------------------------------
426      !
427      vt_i (:,:) = 0._wp
428      DO jl = 1, jpl
429         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
430      END DO
431
432      bv_i(:,:) = 0._wp
433      DO jl = 1, jpl
434         DO jk = 1, nlay_i
435            DO jj = 1, jpj
436               DO ji = 1, jpi
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 )
442               END DO
443            END DO
444         END DO
445      END DO
446      !
447   END SUBROUTINE lim_var_bv
448
449
450   SUBROUTINE lim_var_salprof1d( kideb, kiut )
451      !!-------------------------------------------------------------------
452      !!                  ***  ROUTINE lim_thd_salprof1d  ***
453      !!
454      !! ** Purpose :   1d computation of the sea ice salinity profile
455      !!                Works with 1d vectors and is used by thermodynamic modules
456      !!-------------------------------------------------------------------
457      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
458      !
459      INTEGER  ::   ji, jk    ! dummy loop indices
460      INTEGER  ::   ii, ij    ! local integers
461      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars
462      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      -
463      !
464      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
465      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
466      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
467      !!---------------------------------------------------------------------
468
469      CALL wrk_alloc( jpij, z_slope_s )
470
471      !---------------------------------------
472      ! Vertically constant, constant in time
473      !---------------------------------------
474      IF( nn_icesal == 1 )   s_i_1d(:,:) = rn_icesal
475
476      !------------------------------------------------------
477      ! Vertically varying salinity profile, varying in time
478      !------------------------------------------------------
479
480      IF(  nn_icesal == 2  ) THEN
481         !
482         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
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) )
485         END DO
486
487         ! Weighting factor between zs_zero and zs_inf
488         !---------------------------------------------
489         zfac0 = 1._wp / ( zsi0 - zsi1 )
490         zfac1 = zsi1 / ( zsi1 - zsi0 )
491         DO jk = 1, nlay_i
492            DO ji = kideb, kiut
493               ii =  MOD( npb(ji) - 1 , jpi ) + 1
494               ij =     ( npb(ji) - 1 ) / jpi + 1
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
500               ! this is to force a constant salinity profile in the Baltic Sea
501               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
502               !
503               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch )
504               !
505               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
506               ! weighting the profile
507               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
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
512
513      ENDIF 
514
515      !-------------------------------------------------------
516      ! Vertically varying salinity profile, constant in time
517      !-------------------------------------------------------
518
519      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
520         !
521         sm_i_1d(:) = 2.30_wp
522         !
523         DO jk = 1, nlay_i
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 ) ) ) )
526            DO ji = kideb, kiut
527               s_i_1d(ji,jk) = zsal
528            END DO
529         END DO
530         !
531      ENDIF
532      !
533      CALL wrk_dealloc( jpij, z_slope_s )
534      !
535   END SUBROUTINE lim_var_salprof1d
536
537   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
785#else
786   !!----------------------------------------------------------------------
787   !!   Default option         Dummy module          NO  LIM3 sea-ice model
788   !!----------------------------------------------------------------------
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
799   END SUBROUTINE lim_var_bv
800   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
801   END SUBROUTINE lim_var_salprof1d
802   SUBROUTINE lim_var_zapsmall
803   END SUBROUTINE lim_var_zapsmall
804   SUBROUTINE lim_var_itd
805   END SUBROUTINE lim_var_itd
806#endif
807
808   !!======================================================================
809END MODULE limvar
Note: See TracBrowser for help on using the repository browser.