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.
icevar.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90 @ 8496

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

debug icevar

File size: 34.2 KB
Line 
1MODULE icevar
2   !!======================================================================
3   !!                       ***  MODULE icevar ***
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   !!                        - tm_i (jpi,jpj) !mean ice temperature
30   !!======================================================================
31   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
32   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
33   !!            3.5  ! 2012    (M. Vancoppenolle)  add ice_var_itd
34   !!            3.6  ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd
35   !!----------------------------------------------------------------------
36#if defined key_lim3
37   !!----------------------------------------------------------------------
38   !!   'key_lim3'                                      LIM3 sea-ice model
39   !!----------------------------------------------------------------------
40   !!   ice_var_agg       : integrate variables over layers and categories
41   !!   ice_var_glo2eqv   : transform from VGLO to VEQV
42   !!   ice_var_eqv2glo   : transform from VEQV to VGLO
43   !!   ice_var_salprof   : salinity profile in the ice
44   !!   ice_var_bv        : brine volume
45   !!   ice_var_salprof1d : salinity profile in the ice 1D
46   !!   ice_var_zapsmall  : remove very small area and volume
47   !!   ice_var_itd       : convert 1-cat to multiple cat
48   !!----------------------------------------------------------------------
49   USE par_oce        ! ocean parameters
50   USE phycst         ! physical constants (ocean directory)
51   USE sbc_oce , ONLY : sss_m
52   USE ice            ! ice variables
53   USE ice1D          ! ice variables (thermodynamics)
54   !
55   USE in_out_manager ! I/O manager
56   USE lib_mpp        ! MPP library
57   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
58
59   IMPLICIT NONE
60   PRIVATE
61
62   PUBLIC   ice_var_agg         
63   PUBLIC   ice_var_glo2eqv     
64   PUBLIC   ice_var_eqv2glo     
65   PUBLIC   ice_var_salprof     
66   PUBLIC   ice_var_bv           
67   PUBLIC   ice_var_salprof1d   
68   PUBLIC   ice_var_zapsmall
69   PUBLIC   ice_var_itd
70
71   !!----------------------------------------------------------------------
72   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
73   !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
74   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
75   !!----------------------------------------------------------------------
76CONTAINS
77
78   SUBROUTINE ice_var_agg( kn )
79      !!------------------------------------------------------------------
80      !!                ***  ROUTINE ice_var_agg  ***
81      !!
82      !! ** Purpose :   aggregates ice-thickness-category variables to
83      !!              all-ice variables, i.e. it turns VGLO into VAGG
84      !!------------------------------------------------------------------
85      INTEGER, INTENT( in ) ::   kn     ! =1 state variables only
86      !                                 ! >1 state variables + others
87      !
88      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
89      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z1_at_i, z1_vt_i
90      !!------------------------------------------------------------------
91      !
92      !                                      ! integrated values
93      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
94      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
95      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
96      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
97      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
98
99      ! MV MP 2016
100      IF ( ln_pnd ) THEN                     ! Melt pond
101         at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 )
102         vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
103      ENDIF
104      ! END MP 2016
105
106      DO jj = 1, jpj                         ! open water fraction
107         DO ji = 1, jpi
108            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   
109         END DO
110      END DO
111!!gm I think this should do the work :
112!      ato_i(:,:) = MAX( 1._wp - at_i(:,:), 0._wp ) 
113!!gm end
114
115      IF( kn > 1 ) THEN
116         !
117         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) )
118         WHERE( at_i(:,:) > epsi10 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
119         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
120         END WHERE
121         WHERE( vt_i(:,:) > epsi10 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
122         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
123         END WHERE
124         !
125         !                          ! mean ice/snow thickness
126         htm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
127         htm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
128         !         
129         !                          ! mean temperature (K), salinity and age
130         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
131         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
132         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
133         !
134         tm_i (:,:) = 0._wp
135         smt_i(:,:) = 0._wp
136         DO jl = 1, jpl
137            DO jk = 1, nlay_i
138               tm_i (:,:) = tm_i (:,:) + r1_nlay_i * t_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
139               smt_i(:,:) = smt_i(:,:) + r1_nlay_i * s_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
140            END DO
141         END DO
142         !
143!!gm  QUESTION 1 : why salinity is named smt_i  and not just sm_i ?   since the 4D field is named s_i. (NB for temp: tm_i, t_i)
144         !
145         DEALLOCATE( z1_at_i , z1_vt_i )
146      ENDIF
147      !
148   END SUBROUTINE ice_var_agg
149
150
151   SUBROUTINE ice_var_glo2eqv
152      !!------------------------------------------------------------------
153      !!                ***  ROUTINE ice_var_glo2eqv ***
154      !!
155      !! ** Purpose :   computes equivalent variables as function of 
156      !!              global variables, i.e. it turns VGLO into VEQV
157      !!------------------------------------------------------------------
158      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
159      REAL(wp) ::   ze_i, z1_CpR, zdiscrim, zaaa, z1_2aaa             ! local scalars
160      REAL(wp) ::   ze_s, zL_Cp , ztmelts , zbbb, zccc                !   -      -
161      REAL(wp) ::   zhmax, z1_zhmax, zsm_i, zcpMcpic, zt_i   !   -      -
162      REAL(wp) ::   zlay_i, zlay_s   !   -      -
163      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
164      !!------------------------------------------------------------------
165
166!!gm Question 1:  here use epsi20 , in ice_var_agg it is epsi10 which is used...  why ??
167
168!!gm Question 2:  It is possible to define existence of sea-ice in a common way between
169!!                ice area and ice volume ?
170!!                the idea is to be able to define one for all at the begining of this routine
171!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
172
173      !-------------------------------------------------------
174      ! Ice thickness, snow thickness, ice salinity, ice age
175      !-------------------------------------------------------
176      !                                            !--- inverse of the ice area
177      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
178      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
179      END WHERE
180      !
181      ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness
182
183     zhmax    =          hi_max(jpl)
184      z1_zhmax =  1._wp / hi_max(jpl)               
185      WHERE( ht_i(:,:,jpl) > zhmax )               !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area
186         ht_i  (:,:,jpl) = zhmax
187         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
188         z1_a_i(:,:,jpl) = zhmax / v_i(:,:,jpl)          ! NB: v_i always /=0 as ht_i > hi_max
189      END WHERE
190
191      ht_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)    !--- snow thickness
192     
193      o_i(:,:,:)  = oa_i(:,:,:) * z1_a_i(:,:,:)    !--- ice age
194
195      IF( nn_icesal == 2 )THEN                     !--- salinity (with a minimum value imposed everywhere)
196         WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) )
197         ELSEWHERE                      ;   sm_i(:,:,:) = rn_simin
198         END WHERE
199      ENDIF
200
201      CALL ice_var_salprof      ! salinity profile
202
203      !-------------------
204      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere)
205      !-------------------
206      zlay_i   = REAL( nlay_i , wp )    ! number of layers
207      zaaa     = cpic                   ! Conversion q(S,T) -> T (second order equation)
208      z1_2aaa  = 1._wp / ( 2._wp * zaaa )
209      zcpMcpic = rcp - cpic
210      DO jl = 1, jpl
211         DO jk = 1, nlay_i
212            DO jj = 1, jpj
213               DO ji = 1, jpi
214                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area
215                     !
216                     ze_i    =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
217                     ztmelts = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
218                     !
219                     zbbb     =  zcpMcpic * ztmelts + ze_i * r1_rhoic - lfus
220                     zccc     =  lfus * ztmelts
221                     zdiscrim =  SQRT(  MAX( zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp)  )
222                     zt_i     = - ( zbbb + zdiscrim ) * z1_2aaa                          ! [C]
223                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( zt_i , ztmelts )  ) + rt0   ! [K] with bounds: -100 < zt_i < ztmelts
224                     !
225                  ELSE                                !--- no ice
226                     t_i(ji,jj,jk,jl) =  rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C)
227!!clem: I think we should impose rt0 instead
228                  ENDIF
229               END DO
230            END DO
231         END DO
232      END DO
233
234      !--------------------
235      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere)
236      !--------------------
237      z1_CpR = 1._wp / ( cpic * rhosn )
238      zL_Cp  = lfus  /   cpic
239      zlay_s = REAL( nlay_s , wp )
240      DO jk = 1, nlay_s
241         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
242            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0       
243         ELSEWHERE                           !--- no ice
244            t_s(:,:,jk,:) =  rt0 - 100._wp         ! impose 173.15 K (i.e. -100 C)
245         END WHERE
246      END DO
247!!gm perhaps better like this (?) :  (if this compile, yes! one test instead of nlay_s tests)
248!      WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
249!         DO jk = 1, nlay_s
250!            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0       
251!         END DO
252!      ELSEWHERE                           !--- no ice
253!         DO jk = 1, nlay_s
254!            t_s(:,:,jk,:) =  rt0 - 100._wp      ! impose 173.15 K (i.e. -100 C)
255!         END DO
256!      END WHERE
257!!gm end better ?
258
259      ! integrated values
260      vt_i (:,:) = SUM( v_i, dim=3 )
261      vt_s (:,:) = SUM( v_s, dim=3 )
262      at_i (:,:) = SUM( a_i, dim=3 )
263
264! MV MP 2016
265! probably should resum for melt ponds ???
266
267      !
268   END SUBROUTINE ice_var_glo2eqv
269
270
271   SUBROUTINE ice_var_eqv2glo
272      !!------------------------------------------------------------------
273      !!                ***  ROUTINE ice_var_eqv2glo ***
274      !!
275      !! ** Purpose :   computes global variables as function of
276      !!              equivalent variables,  i.e. it turns VEQV into VGLO
277      !!------------------------------------------------------------------
278      !
279      v_i  (:,:,:) = ht_i(:,:,:) * a_i(:,:,:)
280      v_s  (:,:,:) = ht_s(:,:,:) * a_i(:,:,:)
281      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
282      !
283   END SUBROUTINE ice_var_eqv2glo
284
285
286   SUBROUTINE ice_var_salprof
287      !!------------------------------------------------------------------
288      !!                ***  ROUTINE ice_var_salprof ***
289      !!
290      !! ** Purpose :   computes salinity profile in function of bulk salinity     
291      !!
292      !! ** Method  : If bulk salinity greater than zsi1,
293      !!              the profile is assumed to be constant (S_inf)
294      !!              If bulk salinity lower than zsi0,
295      !!              the profile is linear with 0 at the surface (S_zero)
296      !!              If it is between zsi0 and zsi1, it is a
297      !!              alpha-weighted linear combination of s_inf and s_zero
298      !!
299      !! ** References : Vancoppenolle et al., 2007
300      !!------------------------------------------------------------------
301      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
302      REAL(wp) ::   zsal, z1_dS
303      REAL(wp) ::   zargtemp , zs0, zs
304      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
305      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
306      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
307      !!------------------------------------------------------------------
308
309!!gm  much much more secure to defined when reading nn_icesal in the namelist integers =1, 2, 3  with explicit names
310!!       for example np_Scst_noProfile = 1  ;  np_Svar_linProfile = 2   ;   np_Scst_fixProfile
311
312!!gm Question: Remove the option 3 ?  How many years since it last use ?
313
314      SELECT CASE ( nn_icesal )
315      !
316      !              !---------------------------------------!
317      CASE( 1 )      !  constant salinity in time and space  !
318         !           !---------------------------------------!
319         s_i (:,:,:,:) = rn_icesal
320         sm_i(:,:,:)   = rn_icesal
321         !
322         !           !---------------------------------------------!
323      CASE( 2 )      !  time varying salinity with linear profile  !
324         !           !---------------------------------------------!
325         !
326         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
327         !
328         DO jk = 1, nlay_i
329            s_i(:,:,jk,:)  = sm_i(:,:,:)
330         END DO
331         !                                      ! Slope of the linear profile
332         WHERE( ht_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * sm_i(:,:,:) / ht_i(:,:,:)
333         ELSEWHERE                       ;   z_slope_s(:,:,:) = 0._wp
334         END WHERE
335         !
336         z1_dS = 1._wp / ( zsi1 - zsi0 )
337         DO jl = 1, jpl
338            DO jj = 1, jpj
339               DO ji = 1, jpi
340                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - sm_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
341                  !                             ! force a constant profile when SSS too low (Baltic Sea)
342                  IF( 2._wp * sm_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
343               END DO
344            END DO
345         END DO
346
347         ! Computation of the profile
348         DO jl = 1, jpl
349            DO jk = 1, nlay_i
350               DO jj = 1, jpj
351                  DO ji = 1, jpi
352                     !                          ! linear profile with 0 surface value
353                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
354                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)     ! weighting the profile
355                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
356                  END DO
357               END DO
358            END DO
359         END DO
360         !
361         DEALLOCATE( z_slope_s , zalpha )
362         !
363         !           !-------------------------------------------!
364      CASE( 3 )      ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
365         !           !-------------------------------------------!                                   (mean = 2.30)
366         !
367         sm_i(:,:,:) = 2.30_wp
368!!gm Remark: if we keep the case 3, then compute an store one for all time-step
369!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
370!         DO jk = 1, nlay_i
371!            s_i(:,:,jk,:) = S_prof(jk)
372!         END DO
373!!gm end
374         !
375         DO jl = 1, jpl
376            DO jk = 1, nlay_i
377               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
378               s_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
379            END DO
380         END DO
381         !
382      END SELECT
383      !
384   END SUBROUTINE ice_var_salprof
385
386
387   SUBROUTINE ice_var_bv
388      !!------------------------------------------------------------------
389      !!                ***  ROUTINE ice_var_bv ***
390      !!
391      !! ** Purpose :   computes mean brine volume (%) in sea ice
392      !!
393      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
394      !!
395      !! References : Vancoppenolle et al., JGR, 2007
396      !!------------------------------------------------------------------
397      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
398      !!------------------------------------------------------------------
399      !
400!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
401!!   instead of setting everything to zero as just below
402      bv_i (:,:,:) = 0._wp
403      DO jl = 1, jpl
404         DO jk = 1, nlay_i
405            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
406               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
407            END WHERE
408         END DO
409      END DO
410      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
411      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
412      END WHERE
413      !
414   END SUBROUTINE ice_var_bv
415
416
417   SUBROUTINE ice_var_salprof1d
418      !!-------------------------------------------------------------------
419      !!                  ***  ROUTINE ice_var_salprof1d  ***
420      !!
421      !! ** Purpose :   1d computation of the sea ice salinity profile
422      !!                Works with 1d vectors and is used by thermodynamic modules
423      !!-------------------------------------------------------------------
424      INTEGER  ::   ji, jk    ! dummy loop indices
425      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
426      REAL(wp) ::   zalpha, zs, zs0              !   -      -
427      !
428      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s   !
429      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
430      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
431      !!---------------------------------------------------------------------
432      !
433      SELECT CASE ( nn_icesal )
434      !
435      !              !---------------------------------------!
436      CASE( 1 )      !  constant salinity in time and space  !
437         !           !---------------------------------------!
438         s_i_1d(:,:) = rn_icesal
439         !
440         !           !---------------------------------------------!
441      CASE( 2 )      !  time varying salinity with linear profile  !
442         !           !---------------------------------------------!
443         !
444         ALLOCATE( z_slope_s(jpij) )
445         !
446         DO ji = 1, nidx          ! Slope of the linear profile
447            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
448            z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
449         END DO
450
451         z1_dS = 1._wp / ( zsi1 - zsi0 )
452         DO jk = 1, nlay_i
453            DO ji = 1, nidx
454               zalpha = MAX(  0._wp , MIN(  ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp  )  )
455               IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) )   zalpha = 0._wp               ! cst profile when SSS too low (Baltic Sea)
456               !
457               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i   ! linear profile with 0 surface value
458               zs  = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji)                      ! weighting the profile
459               s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )                     ! bounding salinity
460            END DO
461         END DO
462         !
463         DEALLOCATE( z_slope_s )
464
465         !           !-------------------------------------------!
466      CASE( 3 )      ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
467         !           !-------------------------------------------!                                   (mean = 2.30)
468         !
469         sm_i_1d(:) = 2.30_wp
470         !
471!!gm cf remark in ice_var_salprof routine, CASE( 3 )
472         DO jk = 1, nlay_i
473            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
474            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
475            DO ji = 1, nidx
476               s_i_1d(ji,jk) = zsal
477            END DO
478         END DO
479         !
480      END SELECT
481      !
482   END SUBROUTINE ice_var_salprof1d
483
484
485   SUBROUTINE ice_var_zapsmall
486      !!-------------------------------------------------------------------
487      !!                   ***  ROUTINE ice_var_zapsmall ***
488      !!
489      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
490      !!-------------------------------------------------------------------
491      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
492      REAL(wp) ::   zsal, zvi, zvs, zei, zes, zvp
493      !!-------------------------------------------------------------------
494      !
495      DO jl = 1, jpl       !==  loop over the categories  ==!
496         !
497         !-----------------------------------------------------------------
498         ! Zap ice energy and use ocean heat to melt ice
499         !-----------------------------------------------------------------
500         DO jk = 1, nlay_i
501            DO jj = 1 , jpj
502               DO ji = 1 , jpi
503!!gm I think we can do better/faster  :  to be discussed
504                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
505                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
506                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
507                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
508                  zei              = e_i(ji,jj,jk,jl)
509                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
510                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
511                  ! update exchanges with ocean
512                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
513               END DO
514            END DO
515         END DO
516
517         DO jj = 1 , jpj
518            DO ji = 1 , jpi
519               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
520               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
521               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
522                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
523               zsal = smv_i(ji,jj,  jl)
524               zvi  = v_i  (ji,jj,  jl)
525               zvs  = v_s  (ji,jj,  jl)
526               zes  = e_s  (ji,jj,1,jl)
527               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw )   zvp  = v_ip (ji,jj  ,jl)
528               !-----------------------------------------------------------------
529               ! Zap snow energy
530               !-----------------------------------------------------------------
531               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
532               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
533
534               !-----------------------------------------------------------------
535               ! zap ice and snow volume, add water and salt to ocean
536               !-----------------------------------------------------------------
537               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
538               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
539               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
540               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
541               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
542               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
543               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
544
545               ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch
546               ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch
547
548               ! MV MP 2016
549               IF ( ln_pnd ) THEN
550                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
551                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
552                  IF ( ln_pnd_fw )   wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice
553               ENDIF
554               ! END MV MP 2016
555
556               ! update exchanges with ocean
557               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
558               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
559               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
560               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
561            END DO
562         END DO
563      END DO 
564
565      ! to be sure that at_i is the sum of a_i(jl)
566      at_i (:,:) = a_i(:,:,1)
567     vt_i (:,:) = v_i(:,:,1)
568      DO jl = 2, jpl
569         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
570        vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
571      END DO
572
573      ! open water = 1 if at_i=0 (no re-calculation of ato_i here)
574      DO jj = 1, jpj
575         DO ji = 1, jpi
576            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
577            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
578         END DO
579      END DO
580      !
581   END SUBROUTINE ice_var_zapsmall
582
583
584   SUBROUTINE ice_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
585      !!------------------------------------------------------------------
586      !!                ***  ROUTINE ice_var_itd   ***
587      !!
588      !! ** Purpose :  converting 1-cat ice to multiple ice categories
589      !!
590      !!                  ice thickness distribution follows a gaussian law
591      !!               around the concentration of the most likely ice thickness
592      !!                           (similar as iceist.F90)
593      !!
594      !! ** Method:   Iterative procedure
595      !!               
596      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
597      !!
598      !!               2) Check whether the distribution conserves area and volume, positivity and
599      !!                  category boundaries
600      !!             
601      !!               3) If not (input ice is too thin), the last category is empty and
602      !!                  the number of categories is reduced (jpl-1)
603      !!
604      !!               4) Iterate until ok (SUM(itest(:) = 4)
605      !!
606      !! ** Arguments : zhti: 1-cat ice thickness
607      !!                zhts: 1-cat snow depth
608      !!                zai : 1-cat ice concentration
609      !!
610      !! ** Output    : jpl-cat
611      !!
612      !!  (Example of application: BDY forcings when input are cell averaged) 
613      !!-------------------------------------------------------------------
614      INTEGER  :: ji, jk, jl             ! dummy loop indices
615      INTEGER  :: ijpij, i_fill, jl0 
616      REAL(wp) :: zarg, zV, zconv, zdh, zdv
617      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
618      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
619      INTEGER , DIMENSION(4)                  ::   itest
620      !!-------------------------------------------------------------------
621
622      !--------------------------------------------------------------------
623      ! initialisation of variables
624      !--------------------------------------------------------------------
625      ijpij = SIZE( zhti , 1 )
626      zht_i(1:ijpij,1:jpl) = 0._wp
627      zht_s(1:ijpij,1:jpl) = 0._wp
628      za_i (1:ijpij,1:jpl) = 0._wp
629
630      ! ----------------------------------------
631      ! distribution over the jpl ice categories
632      ! ----------------------------------------
633      DO ji = 1, ijpij
634         
635         IF( zhti(ji) > 0._wp ) THEN
636
637            ! find which category (jl0) the input ice thickness falls into
638            jl0 = jpl
639            DO jl = 1, jpl
640               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
641                  jl0 = jl
642                  CYCLE
643               ENDIF
644            END DO
645
646            ! initialisation of tests
647            itest(:)  = 0
648         
649            i_fill = jpl + 1                                             !====================================
650            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
651               ! iteration                                               !====================================
652               i_fill = i_fill - 1
653               
654               ! initialisation of ice variables for each try
655               zht_i(ji,1:jpl) = 0._wp
656               za_i (ji,1:jpl) = 0._wp
657               itest(:)        = 0     
658               
659               ! *** case very thin ice: fill only category 1
660               IF ( i_fill == 1 ) THEN
661                  zht_i(ji,1) = zhti(ji)
662                  za_i (ji,1) = zai (ji)
663                 
664               ! *** case ice is thicker: fill categories >1
665               ELSE
666
667                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean
668                  DO jl = 1, i_fill - 1
669                     zht_i(ji,jl) = hi_mean(jl)
670                  END DO
671                 
672                  ! Concentrations in the (i_fill-1) categories
673                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
674                  DO jl = 1, i_fill - 1
675                     IF ( jl /= jl0 ) THEN
676                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
677                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
678                     ENDIF
679                  END DO
680                 
681                  ! Concentration in the last (i_fill) category
682                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
683                 
684                  ! Ice thickness in the last (i_fill) category
685                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
686                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
687                 
688                  ! clem: correction if concentration of upper cat is greater than lower cat
689                  !       (it should be a gaussian around jl0 but sometimes it is not)
690                  IF ( jl0 /= jpl ) THEN
691                     DO jl = jpl, jl0+1, -1
692                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
693                           zdv = zht_i(ji,jl) * za_i(ji,jl)
694                           zht_i(ji,jl    ) = 0._wp
695                           za_i (ji,jl    ) = 0._wp
696                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
697                        END IF
698                     ENDDO
699                  ENDIF
700               
701               ENDIF ! case ice is thick or thin
702           
703               !---------------------
704               ! Compatibility tests
705               !---------------------
706               ! Test 1: area conservation
707               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
708               IF ( zconv < epsi06 ) itest(1) = 1
709           
710               ! Test 2: volume conservation
711               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
712               IF ( zconv < epsi06 ) itest(2) = 1
713               
714               ! Test 3: thickness of the last category is in-bounds ?
715               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
716               
717               ! Test 4: positivity of ice concentrations
718               itest(4) = 1
719               DO jl = 1, i_fill
720                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
721               END DO
722               !                                         !============================
723            END DO                                       ! end iteration on categories
724               !                                         !============================
725         ENDIF ! if zhti > 0
726      END DO ! i loop
727
728      ! ------------------------------------------------
729      ! Adding Snow in each category where za_i is not 0
730      ! ------------------------------------------------
731      DO jl = 1, jpl
732         DO ji = 1, ijpij
733            IF( za_i(ji,jl) > 0._wp ) THEN
734               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
735               ! In case snow load is in excess that would lead to transformation from snow to ice
736               ! Then, transfer the snow excess into the ice (different from icethd_dh)
737               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
738               ! recompute ht_i, ht_s avoiding out of bounds values
739               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
740               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
741            ENDIF
742         END DO
743      END DO
744      !
745    END SUBROUTINE ice_var_itd
746
747#else
748   !!----------------------------------------------------------------------
749   !!   Default option         Dummy module          NO  LIM3 sea-ice model
750   !!----------------------------------------------------------------------
751#endif
752
753   !!======================================================================
754END MODULE icevar
Note: See TracBrowser for help on using the repository browser.