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

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

changes in style - part2 -

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