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

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 @ 5059

Last change on this file since 5059 was 5055, checked in by clem, 9 years ago

LIM3: removing par_ice and put jpl, nlay_i and nlay_s as namelist parameters
  • Property svn:keywords set to Id
File size: 32.9 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) - epsi10 ) ) 
127                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch   ! ice salinity
128                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
129                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * 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 = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !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) , epsi10 ) * rswitch
165               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
166               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
167            END DO
168         END DO
169      END DO
170
171      IF(  num_sal == 2  )THEN
172         DO jl = 1, jpl
173            DO jj = 1, jpj
174               DO ji = 1, jpi
175                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !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) , epsi10 ) * rswitch
177               END DO
178            END DO
179         END DO
180      ENDIF
181
182      CALL lim_var_salprof      ! salinity profile
183
184      !-------------------
185      ! Ice temperatures
186      !-------------------
187      DO jl = 1, jpl
188         DO jk = 1, nlay_i
189            DO jj = 1, jpj
190               DO ji = 1, jpi
191                  !                                                              ! Energy of melting q(S,T) [J.m-3]
192                  rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes
193                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 
194                  zq_i    = zq_i * unit_fac                             ! convert units
195                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt              ! Ice layer melt temperature
196                  !
197                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
198                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i / rhoic - lfus
199                  zccc       =  lfus * (ztmelts-rtt)
200                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
201                  t_i(ji,jj,jk,jl) = rtt + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )
202                  t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rtt < t_i < rtt
203               END DO
204            END DO
205         END DO
206      END DO
207
208      !--------------------
209      ! Snow temperatures
210      !--------------------
211      zfac1 = 1._wp / ( rhosn * cpic )
212      zfac2 = lfus / cpic 
213      DO jl = 1, jpl
214         DO jk = 1, nlay_s
215            DO jj = 1, jpj
216               DO ji = 1, jpi
217                  !Energy of melting q(S,T) [J.m-3]
218                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes
219                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)
220                  zq_s  = zq_s * unit_fac                                    ! convert units
221                  !
222                  t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 )
223                  t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rtt < t_i < rtt
224               END DO
225            END DO
226         END DO
227      END DO
228
229      !-------------------
230      ! Mean temperature
231      !-------------------
232      tm_i(:,:) = 0._wp
233      DO jl = 1, jpl
234         DO jk = 1, nlay_i
235            DO jj = 1, jpj
236               DO ji = 1, jpi
237                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
238                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
239                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
240               END DO
241            END DO
242         END DO
243      END DO
244      !
245   END SUBROUTINE lim_var_glo2eqv
246
247
248   SUBROUTINE lim_var_eqv2glo
249      !!------------------------------------------------------------------
250      !!                ***  ROUTINE lim_var_eqv2glo ***
251      !!
252      !! ** Purpose :   computes global variables as function of equivalent variables
253      !!                i.e. it turns VEQV into VGLO
254      !! ** Method  :
255      !!
256      !! ** History :  (01-2006) Martin Vancoppenolle, UCL-ASTR
257      !!------------------------------------------------------------------
258      !
259      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
260      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
261      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
262      oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:)
263      !
264   END SUBROUTINE lim_var_eqv2glo
265
266
267   SUBROUTINE lim_var_salprof
268      !!------------------------------------------------------------------
269      !!                ***  ROUTINE lim_var_salprof ***
270      !!
271      !! ** Purpose :   computes salinity profile in function of bulk salinity     
272      !!
273      !! ** Method  : If bulk salinity greater than zsi1,
274      !!              the profile is assumed to be constant (S_inf)
275      !!              If bulk salinity lower than zsi0,
276      !!              the profile is linear with 0 at the surface (S_zero)
277      !!              If it is between zsi0 and zsi1, it is a
278      !!              alpha-weighted linear combination of s_inf and s_zero
279      !!
280      !! ** References : Vancoppenolle et al., 2007
281      !!------------------------------------------------------------------
282      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
283      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac, zsal
284      REAL(wp) ::   zswi0, zswi01, zswibal, zargtemp , zs_zero   
285      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha
286      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
287      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
288      !!------------------------------------------------------------------
289
290      CALL wrk_alloc( jpi, jpj, jpl, z_slope_s, zalpha )
291
292      !---------------------------------------
293      ! Vertically constant, constant in time
294      !---------------------------------------
295      IF(  num_sal == 1  )   s_i(:,:,:,:) = bulk_sal
296
297      !-----------------------------------
298      ! Salinity profile, varying in time
299      !-----------------------------------
300      IF(  num_sal == 2  ) THEN
301         !
302         DO jk = 1, nlay_i
303            s_i(:,:,jk,:)  = sm_i(:,:,:)
304         END DO
305         !
306         DO jl = 1, jpl                               ! Slope of the linear profile
307            DO jj = 1, jpj
308               DO ji = 1, jpi
309                  z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) )
310               END DO
311            END DO
312         END DO
313         !
314         dummy_fac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
315         dummy_fac1 = zsi1 / ( zsi1 - zsi0 )
316         !
317         zalpha(:,:,:) = 0._wp
318         DO jl = 1, jpl
319            DO jj = 1, jpj
320               DO ji = 1, jpi
321                  ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
322                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , zsi0 - sm_i(ji,jj,jl) ) ) 
323                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
324                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) ) 
325                  ! If 2.sm_i GE sss_m then zswibal = 1
326                  ! this is to force a constant salinity profile in the Baltic Sea
327                  zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
328                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )
329                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal )
330               END DO
331            END DO
332         END DO
333
334         dummy_fac = 1._wp / REAL( nlay_i )                   ! Computation of the profile
335         DO jl = 1, jpl
336            DO jk = 1, nlay_i
337               DO jj = 1, jpj
338                  DO ji = 1, jpi
339                     !                                      ! linear profile with 0 at the surface
340                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac
341                     !                                      ! weighting the profile
342                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
343                  END DO ! ji
344               END DO ! jj
345            END DO ! jk
346         END DO ! jl
347         !
348      ENDIF ! num_sal
349
350      !-------------------------------------------------------
351      ! Vertically varying salinity profile, constant in time
352      !-------------------------------------------------------
353
354      IF(  num_sal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
355         !
356         sm_i(:,:,:) = 2.30_wp
357         !
358         DO jl = 1, jpl
359            DO jk = 1, nlay_i
360               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
361               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
362               s_i(:,:,jk,jl) =  zsal
363            END DO
364         END DO
365         !
366      ENDIF ! num_sal
367      !
368      CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha )
369      !
370   END SUBROUTINE lim_var_salprof
371
372
373   SUBROUTINE lim_var_icetm
374      !!------------------------------------------------------------------
375      !!                ***  ROUTINE lim_var_icetm ***
376      !!
377      !! ** Purpose :   computes mean sea ice temperature
378      !!------------------------------------------------------------------
379      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
380      !!------------------------------------------------------------------
381
382      ! Mean sea ice temperature
383      tm_i(:,:) = 0._wp
384      DO jl = 1, jpl
385         DO jk = 1, nlay_i
386            DO jj = 1, jpj
387               DO ji = 1, jpi
388                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
389                  tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
390                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
391               END DO
392            END DO
393         END DO
394      END DO
395
396   END SUBROUTINE lim_var_icetm
397
398
399   SUBROUTINE lim_var_bv
400      !!------------------------------------------------------------------
401      !!                ***  ROUTINE lim_var_bv ***
402      !!
403      !! ** Purpose :   computes mean brine volume (%) in sea ice
404      !!
405      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
406      !!
407      !! References : Vancoppenolle et al., JGR, 2007
408      !!------------------------------------------------------------------
409      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
410      REAL(wp) ::   zbvi             ! local scalars
411      !!------------------------------------------------------------------
412      !
413      bv_i(:,:) = 0._wp
414      DO jl = 1, jpl
415         DO jk = 1, nlay_i
416            DO jj = 1, jpj
417               DO ji = 1, jpi
418                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  )
419                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   &
420                     &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp)
421                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
422                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi10 )
423               END DO
424            END DO
425         END DO
426      END DO
427      !
428   END SUBROUTINE lim_var_bv
429
430
431   SUBROUTINE lim_var_salprof1d( kideb, kiut )
432      !!-------------------------------------------------------------------
433      !!                  ***  ROUTINE lim_thd_salprof1d  ***
434      !!
435      !! ** Purpose :   1d computation of the sea ice salinity profile
436      !!                Works with 1d vectors and is used by thermodynamic modules
437      !!-------------------------------------------------------------------
438      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
439      !
440      INTEGER  ::   ji, jk    ! dummy loop indices
441      INTEGER  ::   ii, ij    ! local integers
442      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars
443      REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, zs_zero              !   -      -
444      !
445      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
446      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
447      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
448      !!---------------------------------------------------------------------
449
450      CALL wrk_alloc( jpij, z_slope_s )
451
452      !---------------------------------------
453      ! Vertically constant, constant in time
454      !---------------------------------------
455      IF( num_sal == 1 )   s_i_1d(:,:) = bulk_sal
456
457      !------------------------------------------------------
458      ! Vertically varying salinity profile, varying in time
459      !------------------------------------------------------
460
461      IF(  num_sal == 2  ) THEN
462         !
463         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
464            z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) )
465         END DO
466
467         ! Weighting factor between zs_zero and zs_inf
468         !---------------------------------------------
469         dummy_fac0 = 1._wp / ( zsi0 - zsi1 )
470         dummy_fac1 = zsi1 / ( zsi1 - zsi0 )
471         dummy_fac2 = 1._wp / REAL(nlay_i,wp)
472
473         DO jk = 1, nlay_i
474            DO ji = kideb, kiut
475               ii =  MOD( npb(ji) - 1 , jpi ) + 1
476               ij =     ( npb(ji) - 1 ) / jpi + 1
477               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
478               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
479               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws
480               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
481               ! if 2.sm_i GE sss_m then zswibal = 1
482               ! this is to force a constant salinity profile in the Baltic Sea
483               zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
484               !
485               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zswibal )
486               !
487               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2
488               ! weighting the profile
489               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
490            END DO
491         END DO
492
493      ENDIF 
494
495      !-------------------------------------------------------
496      ! Vertically varying salinity profile, constant in time
497      !-------------------------------------------------------
498
499      IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
500         !
501         sm_i_1d(:) = 2.30_wp
502         !
503         DO jk = 1, nlay_i
504            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
505            zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
506            DO ji = kideb, kiut
507               s_i_1d(ji,jk) = zsal
508            END DO
509         END DO
510         !
511      ENDIF
512      !
513      CALL wrk_dealloc( jpij, z_slope_s )
514      !
515   END SUBROUTINE lim_var_salprof1d
516
517   SUBROUTINE lim_var_zapsmall
518      !!-------------------------------------------------------------------
519      !!                   ***  ROUTINE lim_var_zapsmall ***
520      !!
521      !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes
522      !!
523      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
524      !!-------------------------------------------------------------------
525      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
526
527      REAL(wp) ::   zsal, zvi, zvs, zei, zes
528      !!-------------------------------------------------------------------
529
530      DO jl = 1, jpl
531
532         !-----------------------------------------------------------------
533         ! Zap ice energy and use ocean heat to melt ice
534         !-----------------------------------------------------------------
535         DO jk = 1, nlay_i
536            DO jj = 1 , jpj
537               DO ji = 1 , jpi
538                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
539                  zei              = e_i(ji,jj,jk,jl)
540                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
541                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rtt * ( 1._wp - rswitch )
542                  ! update exchanges with ocean
543                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
544               END DO
545            END DO
546         END DO
547
548         DO jj = 1 , jpj
549            DO ji = 1 , jpi
550               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
551               
552               zsal = smv_i(ji,jj,  jl)
553               zvi  = v_i  (ji,jj,  jl)
554               zvs  = v_s  (ji,jj,  jl)
555               zes  = e_s  (ji,jj,1,jl)
556               !-----------------------------------------------------------------
557               ! Zap snow energy
558               !-----------------------------------------------------------------
559               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rtt * ( 1._wp - rswitch )
560               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
561
562               !-----------------------------------------------------------------
563               ! zap ice and snow volume, add water and salt to ocean
564               !-----------------------------------------------------------------
565               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
566               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
567               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
568               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
569               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
570               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
571               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
572
573               ! ice salinity must stay in bounds
574               IF(  num_sal == 2  ) THEN
575                  smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )
576               ENDIF
577               ! update exchanges with ocean
578               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
579               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
580               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
581               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
582            END DO
583         END DO
584      END DO ! jl
585
586      ! to be sure that at_i is the sum of a_i(jl)
587      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
588      !
589   END SUBROUTINE lim_var_zapsmall
590
591   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
592      !!------------------------------------------------------------------
593      !!                ***  ROUTINE lim_var_itd   ***
594      !!
595      !! ** Purpose :  converting 1-cat ice to multiple ice categories
596      !!
597      !!                  ice thickness distribution follows a gaussian law
598      !!               around the concentration of the most likely ice thickness
599      !!                           (similar as limistate.F90)
600      !!
601      !! ** Method:   Iterative procedure
602      !!               
603      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
604      !!
605      !!               2) Check whether the distribution conserves area and volume, positivity and
606      !!                  category boundaries
607      !!             
608      !!               3) If not (input ice is too thin), the last category is empty and
609      !!                  the number of categories is reduced (jpl-1)
610      !!
611      !!               4) Iterate until ok (SUM(itest(:) = 4)
612      !!
613      !! ** Arguments : zhti: 1-cat ice thickness
614      !!                zhts: 1-cat snow depth
615      !!                zai : 1-cat ice concentration
616      !!
617      !! ** Output    : jpl-cat
618      !!
619      !!  (Example of application: BDY forcings when input are cell averaged) 
620      !!
621      !!-------------------------------------------------------------------
622      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
623      !!                    2014    (C. Rousset)        Rewriting
624      !!-------------------------------------------------------------------
625      !! Local variables
626      INTEGER  :: ji, jk, jl             ! dummy loop indices
627      INTEGER  :: ijpij, i_fill, jl0 
628      REAL(wp) :: zarg, zV, zconv, zdh
629      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
630      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
631      INTEGER , POINTER, DIMENSION(:)         ::   itest
632 
633      CALL wrk_alloc( 4, itest )
634      !--------------------------------------------------------------------
635      ! initialisation of variables
636      !--------------------------------------------------------------------
637      ijpij = SIZE(zhti,1)
638      zht_i(1:ijpij,1:jpl) = 0._wp
639      zht_s(1:ijpij,1:jpl) = 0._wp
640      za_i (1:ijpij,1:jpl) = 0._wp
641
642      ! ----------------------------------------
643      ! distribution over the jpl ice categories
644      ! ----------------------------------------
645      DO ji = 1, ijpij
646         
647         IF( zhti(ji) > 0._wp ) THEN
648
649         ! initialisation of tests
650         itest(:)  = 0
651         
652         i_fill = jpl + 1                                             !====================================
653         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
654            ! iteration                                               !====================================
655            i_fill = i_fill - 1
656           
657            ! initialisation of ice variables for each try
658            zht_i(ji,1:jpl) = 0._wp
659            za_i (ji,1:jpl) = 0._wp
660           
661            ! *** case very thin ice: fill only category 1
662            IF ( i_fill == 1 ) THEN
663               zht_i(ji,1) = zhti(ji)
664               za_i (ji,1) = zai (ji)
665
666            ! *** case ice is thicker: fill categories >1
667            ELSE
668
669               ! Fill ice thicknesses except the last one (i_fill) by hmean
670               DO jl = 1, i_fill - 1
671                  zht_i(ji,jl) = hi_mean(jl)
672               END DO
673               
674               ! find which category (jl0) the input ice thickness falls into
675               jl0 = i_fill
676               DO jl = 1, i_fill
677                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
678                     jl0 = jl
679           CYCLE
680                  ENDIF
681               END DO
682               
683               ! Concentrations in the (i_fill-1) categories
684               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
685               DO jl = 1, i_fill - 1
686                  IF ( jl == jl0 ) CYCLE
687                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
688                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
689               END DO
690               
691               ! Concentration in the last (i_fill) category
692               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
693               
694               ! Ice thickness in the last (i_fill) category
695               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
696               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
697               
698            ENDIF ! case ice is thick or thin
699           
700            !---------------------
701            ! Compatibility tests
702            !---------------------
703            ! Test 1: area conservation
704            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
705            IF ( zconv < epsi06 ) itest(1) = 1
706           
707            ! Test 2: volume conservation
708            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
709            IF ( zconv < epsi06 ) itest(2) = 1
710           
711            ! Test 3: thickness of the last category is in-bounds ?
712            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
713           
714            ! Test 4: positivity of ice concentrations
715            itest(4) = 1
716            DO jl = 1, i_fill
717               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
718            END DO           
719                                                           !============================
720         END DO                                            ! end iteration on categories
721                                                           !============================
722         ENDIF ! if zhti > 0
723      END DO ! i loop
724
725      ! ------------------------------------------------
726      ! Adding Snow in each category where za_i is not 0
727      ! ------------------------------------------------
728      DO jl = 1, jpl
729         DO ji = 1, ijpij
730            IF( za_i(ji,jl) > 0._wp ) THEN
731               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
732               ! In case snow load is in excess that would lead to transformation from snow to ice
733               ! Then, transfer the snow excess into the ice (different from limthd_dh)
734               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
735               ! recompute ht_i, ht_s avoiding out of bounds values
736               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
737               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
738            ENDIF
739         ENDDO
740      ENDDO
741
742      CALL wrk_dealloc( 4, itest )
743      !
744    END SUBROUTINE lim_var_itd
745
746
747#else
748   !!----------------------------------------------------------------------
749   !!   Default option         Dummy module          NO  LIM3 sea-ice model
750   !!----------------------------------------------------------------------
751CONTAINS
752   SUBROUTINE lim_var_agg          ! Empty routines
753   END SUBROUTINE lim_var_agg
754   SUBROUTINE lim_var_glo2eqv      ! Empty routines
755   END SUBROUTINE lim_var_glo2eqv
756   SUBROUTINE lim_var_eqv2glo      ! Empty routines
757   END SUBROUTINE lim_var_eqv2glo
758   SUBROUTINE lim_var_salprof      ! Empty routines
759   END SUBROUTINE lim_var_salprof
760   SUBROUTINE lim_var_bv           ! Emtpy routines
761   END SUBROUTINE lim_var_bv
762   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
763   END SUBROUTINE lim_var_salprof1d
764   SUBROUTINE lim_var_zapsmall
765   END SUBROUTINE lim_var_zapsmall
766   SUBROUTINE lim_var_itd
767   END SUBROUTINE lim_var_itd
768#endif
769
770   !!======================================================================
771END MODULE limvar
Note: See TracBrowser for help on using the repository browser.