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

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

changes in style - part3 - move output into proper files and correct a bug in debug mode

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