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

Last change on this file since 8486 was 8486, checked in by clem, 3 years ago

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

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