New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icevar.F90 in NEMO/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icevar.F90 @ 9596

Last change on this file since 9596 was 9570, checked in by nicolasmartin, 6 years ago

Global renaming for core routines (./NEMO)

  • Folders
    • LIM_SRC_3 -> ICE_SRC
    • OPA_SRC -> OCE_SRC
  • CPP key: key_lim3 -> key_si3
  • Modules, (sub)routines and variables names
    • MPI: mpi_comm_opa -> mpi_comm_oce, MPI_COMM_OPA -> MPI_COMM_OCE, mpi_init_opa -> mpi_init_oce
    • AGRIF: agrif_opa_* -> agrif_oce_*, agrif_lim3_* -> agrif_si3_* and few more
    • TOP-PISCES: p.zlim -> p.zice, namp.zlim -> namp.zice
  • Comments
    • NEMO/OPA -> NEMO/OCE
    • ESIM|LIM3 -> SI3
File size: 39.3 KB
Line 
1MODULE icevar
2   !!======================================================================
3   !!                       ***  MODULE icevar ***
4   !!   sea-ice:     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   !!                        - sv_i(jpi,jpj,jpl)
15   !!                        - oa_i(jpi,jpj,jpl)
16   !!                 VEQV : equivalent variables sometimes used in the model
17   !!                        - h_i(jpi,jpj,jpl)
18   !!                        - h_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   !!                        - sm_i(jpi,jpj)  mean ice salinity
29   !!                        - tm_i(jpi,jpj)  mean ice temperature
30   !!                        - tm_s(jpi,jpj)  mean snw temperature
31   !!======================================================================
32   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
33   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
34   !!            3.5  ! 2012    (M. Vancoppenolle)  add ice_var_itd
35   !!            3.6  ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd
36   !!----------------------------------------------------------------------
37#if defined key_si3
38   !!----------------------------------------------------------------------
39   !!   'key_si3'                                       SI3 sea-ice model
40   !!----------------------------------------------------------------------
41   !!   ice_var_agg       : integrate variables over layers and categories
42   !!   ice_var_glo2eqv   : transform from VGLO to VEQV
43   !!   ice_var_eqv2glo   : transform from VEQV to VGLO
44   !!   ice_var_salprof   : salinity profile in the ice
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 jpl-cat
48   !!   ice_var_itd2      : convert N-cat to jpl-cat
49   !!   ice_var_bv        : brine volume
50   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature
51   !!----------------------------------------------------------------------
52   USE dom_oce        ! ocean space and time domain
53   USE phycst         ! physical constants (ocean directory)
54   USE sbc_oce , ONLY : sss_m
55   USE ice            ! sea-ice: variables
56   USE ice1D          ! sea-ice: thermodynamics variables
57   !
58   USE in_out_manager ! I/O manager
59   USE lib_mpp        ! MPP library
60   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
61
62   IMPLICIT NONE
63   PRIVATE
64
65   PUBLIC   ice_var_agg         
66   PUBLIC   ice_var_glo2eqv     
67   PUBLIC   ice_var_eqv2glo     
68   PUBLIC   ice_var_salprof     
69   PUBLIC   ice_var_salprof1d   
70   PUBLIC   ice_var_zapsmall
71   PUBLIC   ice_var_itd
72   PUBLIC   ice_var_itd2
73   PUBLIC   ice_var_bv           
74   PUBLIC   ice_var_enthalpy           
75
76   !!----------------------------------------------------------------------
77   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
78   !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
79   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
80   !!----------------------------------------------------------------------
81CONTAINS
82
83   SUBROUTINE ice_var_agg( kn )
84      !!-------------------------------------------------------------------
85      !!                ***  ROUTINE ice_var_agg  ***
86      !!
87      !! ** Purpose :   aggregates ice-thickness-category variables to
88      !!              all-ice variables, i.e. it turns VGLO into VAGG
89      !!-------------------------------------------------------------------
90      INTEGER, INTENT( in ) ::   kn     ! =1 state variables only
91      !                                 ! >1 state variables + others
92      !
93      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
94      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z1_at_i, z1_vt_i, z1_vt_s
95      !!-------------------------------------------------------------------
96      !
97      !                                      ! integrated values
98      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
99      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
100      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
101      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
102      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
103      !
104      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
105      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
106      !
107      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction 
108
109      IF( kn > 1 ) THEN
110         !
111         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
112         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
113         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
114         END WHERE
115         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
116         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
117         END WHERE
118         WHERE( vt_s(:,:) > epsi20 )   ;   z1_vt_s(:,:) = 1._wp / vt_s(:,:)
119         ELSEWHERE                     ;   z1_vt_s(:,:) = 0._wp
120         END WHERE
121         !
122         !                          ! mean ice/snow thickness
123         hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
124         hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
125         !         
126         !                          ! mean temperature (K), salinity and age
127         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
128         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
129         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
130         sm_i (:,:) = SUM( sv_i(:,:,:)              , dim=3 ) * z1_vt_i(:,:)
131         !
132         tm_i(:,:) = 0._wp
133         tm_s(:,:) = 0._wp
134         DO jl = 1, jpl
135            DO jk = 1, nlay_i
136               tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
137            END DO
138            DO jk = 1, nlay_s
139               tm_s(:,:) = tm_s(:,:) + r1_nlay_s * t_s (:,:,jk,jl) * v_s(:,:,jl) * z1_vt_s(:,:)
140            END DO
141         END DO
142         !
143         !                           ! put rt0 where there is no ice
144         WHERE( at_i(:,:)<=epsi20 )
145            tm_su(:,:) = rt0
146            tm_si(:,:) = rt0
147            tm_i (:,:) = rt0
148            tm_s (:,:) = rt0
149         END WHERE
150
151         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
152      ENDIF
153      !
154   END SUBROUTINE ice_var_agg
155
156
157   SUBROUTINE ice_var_glo2eqv
158      !!-------------------------------------------------------------------
159      !!                ***  ROUTINE ice_var_glo2eqv ***
160      !!
161      !! ** Purpose :   computes equivalent variables as function of 
162      !!              global variables, i.e. it turns VGLO into VEQV
163      !!-------------------------------------------------------------------
164      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
165      REAL(wp) ::   ze_i             ! local scalars
166      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
167      REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
168      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
169      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
170      !!-------------------------------------------------------------------
171
172!!gm Question 2:  It is possible to define existence of sea-ice in a common way between
173!!                ice area and ice volume ?
174!!                the idea is to be able to define one for all at the begining of this routine
175!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
176
177      !---------------------------------------------------------------
178      ! Ice thickness, snow thickness, ice salinity, ice age and ponds
179      !---------------------------------------------------------------
180      !                                            !--- inverse of the ice area
181      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
182      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
183      END WHERE
184      !
185      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
186      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
187      END WHERE
188      !                                           !--- ice thickness
189      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
190
191      zhmax    =          hi_max(jpl)
192      z1_zhmax =  1._wp / hi_max(jpl)               
193      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
194         h_i  (:,:,jpl) = zhmax
195         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
196         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
197      END WHERE
198      !                                           !--- snow thickness
199      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
200      !                                           !--- ice age     
201      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
202      !                                           !--- pond fraction and thickness     
203      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
204      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
205      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
206      END WHERE
207      !
208      !                                           !---  salinity (with a minimum value imposed everywhere)     
209      IF( nn_icesal == 2 ) THEN
210         WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
211         ELSEWHERE                      ;   s_i(:,:,:) = rn_simin
212         END WHERE
213      ENDIF
214      CALL ice_var_salprof   ! salinity profile
215
216      !-------------------
217      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.))
218      !-------------------
219      zlay_i   = REAL( nlay_i , wp )    ! number of layers
220      DO jl = 1, jpl
221         DO jk = 1, nlay_i
222            DO jj = 1, jpj
223               DO ji = 1, jpi
224                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area
225                     !
226                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
227                     ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
228                     ! Conversion q(S,T) -> T (second order equation)
229                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
230                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
231                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
232                     !
233                  ELSE                                !--- no ice
234                     t_i(ji,jj,jk,jl) = rt0
235                  ENDIF
236               END DO
237            END DO
238         END DO
239      END DO
240
241      !--------------------
242      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.))
243      !--------------------
244      zlay_s = REAL( nlay_s , wp )
245      DO jk = 1, nlay_s
246         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
247            t_s(:,:,jk,:) = rt0 + MAX( -100._wp ,  &
248                 &                MIN( r1_cpic * ( -r1_rhosn * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + lfus ) , 0._wp ) )
249         ELSEWHERE                           !--- no ice
250            t_s(:,:,jk,:) = rt0
251         END WHERE
252      END DO
253      !
254      ! integrated values
255      vt_i (:,:) = SUM( v_i, dim=3 )
256      vt_s (:,:) = SUM( v_s, dim=3 )
257      at_i (:,:) = SUM( a_i, dim=3 )
258      !
259   END SUBROUTINE ice_var_glo2eqv
260
261
262   SUBROUTINE ice_var_eqv2glo
263      !!-------------------------------------------------------------------
264      !!                ***  ROUTINE ice_var_eqv2glo ***
265      !!
266      !! ** Purpose :   computes global variables as function of
267      !!              equivalent variables,  i.e. it turns VEQV into VGLO
268      !!-------------------------------------------------------------------
269      !
270      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
271      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
272      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
273      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
274      !
275   END SUBROUTINE ice_var_eqv2glo
276
277
278   SUBROUTINE ice_var_salprof
279      !!-------------------------------------------------------------------
280      !!                ***  ROUTINE ice_var_salprof ***
281      !!
282      !! ** Purpose :   computes salinity profile in function of bulk salinity     
283      !!
284      !! ** Method  : If bulk salinity greater than zsi1,
285      !!              the profile is assumed to be constant (S_inf)
286      !!              If bulk salinity lower than zsi0,
287      !!              the profile is linear with 0 at the surface (S_zero)
288      !!              If it is between zsi0 and zsi1, it is a
289      !!              alpha-weighted linear combination of s_inf and s_zero
290      !!
291      !! ** References : Vancoppenolle et al., 2007
292      !!-------------------------------------------------------------------
293      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
294      REAL(wp) ::   zsal, z1_dS
295      REAL(wp) ::   zargtemp , zs0, zs
296      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
297      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
298      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
299      !!-------------------------------------------------------------------
300
301!!gm Question: Remove the option 3 ?  How many years since it last use ?
302
303      SELECT CASE ( nn_icesal )
304      !
305      !               !---------------------------------------!
306      CASE( 1 )       !  constant salinity in time and space  !
307         !            !---------------------------------------!
308         sz_i(:,:,:,:) = rn_icesal
309         s_i (:,:,:)   = rn_icesal
310         !
311         !            !---------------------------------------------!
312      CASE( 2 )       !  time varying salinity with linear profile  !
313         !            !---------------------------------------------!
314         !
315         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
316         !
317         DO jl = 1, jpl
318            DO jk = 1, nlay_i
319               sz_i(:,:,jk,jl)  = s_i(:,:,jl)
320            END DO
321         END DO
322         !                                      ! Slope of the linear profile
323         WHERE( h_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
324         ELSEWHERE                      ;   z_slope_s(:,:,:) = 0._wp
325         END WHERE
326         !
327         z1_dS = 1._wp / ( zsi1 - zsi0 )
328         DO jl = 1, jpl
329            DO jj = 1, jpj
330               DO ji = 1, jpi
331                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
332                  !                             ! force a constant profile when SSS too low (Baltic Sea)
333                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp 
334               END DO
335            END DO
336         END DO
337         !
338         ! Computation of the profile
339         DO jl = 1, jpl
340            DO jk = 1, nlay_i
341               DO jj = 1, jpj
342                  DO ji = 1, jpi
343                     !                          ! linear profile with 0 surface value
344                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
345                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
346                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
347                  END DO
348               END DO
349            END DO
350         END DO
351         !
352         DEALLOCATE( z_slope_s , zalpha )
353         !
354         !            !-------------------------------------------!
355      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
356         !            !-------------------------------------------!                                   (mean = 2.30)
357         !
358         s_i(:,:,:) = 2.30_wp
359!!gm Remark: if we keep the case 3, then compute an store one for all time-step
360!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
361!         DO jk = 1, nlay_i
362!            sz_i(:,:,jk,:) = S_prof(jk)
363!         END DO
364!!gm end
365         !
366         DO jl = 1, jpl
367            DO jk = 1, nlay_i
368               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
369               sz_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
370            END DO
371         END DO
372         !
373      END SELECT
374      !
375   END SUBROUTINE ice_var_salprof
376
377
378   SUBROUTINE ice_var_salprof1d
379      !!-------------------------------------------------------------------
380      !!                  ***  ROUTINE ice_var_salprof1d  ***
381      !!
382      !! ** Purpose :   1d computation of the sea ice salinity profile
383      !!                Works with 1d vectors and is used by thermodynamic modules
384      !!-------------------------------------------------------------------
385      INTEGER  ::   ji, jk    ! dummy loop indices
386      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
387      REAL(wp) ::   zs, zs0              !   -      -
388      !
389      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
390      REAL(wp), PARAMETER :: zsi0 = 3.5_wp
391      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
392      !!-------------------------------------------------------------------
393      !
394      SELECT CASE ( nn_icesal )
395      !
396      !               !---------------------------------------!
397      CASE( 1 )       !  constant salinity in time and space  !
398         !            !---------------------------------------!
399         sz_i_1d(1:npti,:) = rn_icesal
400         !
401         !            !---------------------------------------------!
402      CASE( 2 )       !  time varying salinity with linear profile  !
403         !            !---------------------------------------------!
404         !
405         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
406         !
407         !                                      ! Slope of the linear profile
408         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
409         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp
410         END WHERE
411         
412         z1_dS = 1._wp / ( zsi1 - zsi0 )
413         DO ji = 1, npti
414            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
415            !                             ! force a constant profile when SSS too low (Baltic Sea)
416            IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
417         END DO
418         !
419         ! Computation of the profile
420         DO jk = 1, nlay_i
421            DO ji = 1, npti
422               !                          ! linear profile with 0 surface value
423               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
424               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
425               sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
426            END DO
427         END DO
428         !
429         DEALLOCATE( z_slope_s, zalpha )
430
431         !            !-------------------------------------------!
432      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
433         !            !-------------------------------------------!                                   (mean = 2.30)
434         !
435         s_i_1d(1:npti) = 2.30_wp
436         !
437!!gm cf remark in ice_var_salprof routine, CASE( 3 )
438         DO jk = 1, nlay_i
439            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
440            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
441            DO ji = 1, npti
442               sz_i_1d(ji,jk) = zsal
443            END DO
444         END DO
445         !
446      END SELECT
447      !
448   END SUBROUTINE ice_var_salprof1d
449
450
451   SUBROUTINE ice_var_zapsmall
452      !!-------------------------------------------------------------------
453      !!                   ***  ROUTINE ice_var_zapsmall ***
454      !!
455      !! ** Purpose :   Remove too small sea ice areas and correct fluxes
456      !!-------------------------------------------------------------------
457      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
458      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
459      !!-------------------------------------------------------------------
460      !
461      DO jl = 1, jpl       !==  loop over the categories  ==!
462         !
463         !-----------------------------------------------------------------
464         ! Zap ice energy and use ocean heat to melt ice
465         !-----------------------------------------------------------------
466         WHERE( a_i(:,:,jl) > epsi10 )   ;   h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
467         ELSEWHERE                       ;   h_i(:,:,jl) = 0._wp
468         END WHERE
469         !
470         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
471         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp
472         END WHERE
473         !
474         DO jk = 1, nlay_i
475            DO jj = 1 , jpj
476               DO ji = 1 , jpi
477                  ! update exchanges with ocean
478                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
479                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
480                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
481               END DO
482            END DO
483         END DO
484         !
485         DO jk = 1, nlay_s
486            DO jj = 1 , jpj
487               DO ji = 1 , jpi
488                  ! update exchanges with ocean
489                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
490                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
491                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
492               END DO
493            END DO
494         END DO
495         !
496         DO jj = 1 , jpj
497            DO ji = 1 , jpi
498               ! update exchanges with ocean
499               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_rdtice
500               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice
501               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
502               !
503               !-----------------------------------------------------------------
504               ! zap ice and snow volume, add water and salt to ocean
505               !-----------------------------------------------------------------
506               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
507               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
508               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
509               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
510               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
511               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
512               !
513               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
514               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
515               !
516               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
517               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
518               !
519            END DO
520         END DO
521         !
522      END DO 
523
524      ! to be sure that at_i is the sum of a_i(jl)
525      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
526      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
527
528      ! open water = 1 if at_i=0
529      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
530      !
531   END SUBROUTINE ice_var_zapsmall
532
533
534   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
535      !!-------------------------------------------------------------------
536      !!                ***  ROUTINE ice_var_itd   ***
537      !!
538      !! ** Purpose :  converting 1-cat ice to multiple ice categories
539      !!
540      !!                  ice thickness distribution follows a gaussian law
541      !!               around the concentration of the most likely ice thickness
542      !!                           (similar as iceistate.F90)
543      !!
544      !! ** Method:   Iterative procedure
545      !!               
546      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
547      !!
548      !!               2) Check whether the distribution conserves area and volume, positivity and
549      !!                  category boundaries
550      !!             
551      !!               3) If not (input ice is too thin), the last category is empty and
552      !!                  the number of categories is reduced (jpl-1)
553      !!
554      !!               4) Iterate until ok (SUM(itest(:) = 4)
555      !!
556      !! ** Arguments : zhti: 1-cat ice thickness
557      !!                zhts: 1-cat snow depth
558      !!                zati: 1-cat ice concentration
559      !!
560      !! ** Output    : jpl-cat
561      !!
562      !!  (Example of application: BDY forcings when input are cell averaged) 
563      !!-------------------------------------------------------------------
564      INTEGER  :: ji, jk, jl             ! dummy loop indices
565      INTEGER  :: idim, i_fill, jl0 
566      REAL(wp) :: zarg, zV, zconv, zdh, zdv
567      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
568      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i ! output ice/snow variables
569      INTEGER , DIMENSION(4)                  ::   itest
570      !!-------------------------------------------------------------------
571      !
572      ! ----------------------------------------
573      ! distribution over the jpl ice categories
574      ! ----------------------------------------
575      ! a gaussian distribution for ice concentration is used
576      ! then we check whether the distribution fullfills
577      ! volume and area conservation, positivity and ice categories bounds
578      idim = SIZE( zhti , 1 )
579      zh_i(1:idim,1:jpl) = 0._wp
580      zh_s(1:idim,1:jpl) = 0._wp
581      za_i(1:idim,1:jpl) = 0._wp
582      !
583      DO ji = 1, idim
584         !
585         IF( zhti(ji) > 0._wp ) THEN
586            !
587            ! find which category (jl0) the input ice thickness falls into
588            jl0 = jpl
589            DO jl = 1, jpl
590               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
591                  jl0 = jl
592                  CYCLE
593               ENDIF
594            END DO
595            !
596            itest(:) = 0
597            i_fill   = jpl + 1                                            !------------------------------------
598            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
599               !                                                          !------------------------------------
600               i_fill = i_fill - 1
601               !
602               zh_i(ji,1:jpl) = 0._wp
603               za_i(ji,1:jpl) = 0._wp
604               itest(:)       = 0     
605               !
606               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
607                  zh_i(ji,1) = zhti(ji)
608                  za_i (ji,1) = zati (ji)
609               ELSE                         !-- case ice is thicker: fill categories >1
610                  ! thickness
611                  DO jl = 1, i_fill - 1
612                     zh_i(ji,jl) = hi_mean(jl)
613                  END DO
614                  !
615                  ! concentration
616                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
617                  DO jl = 1, i_fill - 1
618                     IF ( jl /= jl0 ) THEN
619                        zarg        = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
620                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
621                     ENDIF
622                  END DO
623                  !
624                  ! last category
625                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
626                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
627                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
628                  !
629                  ! correction if concentration of upper cat is greater than lower cat
630                  !    (it should be a gaussian around jl0 but sometimes it is not)
631                  IF ( jl0 /= jpl ) THEN
632                     DO jl = jpl, jl0+1, -1
633                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
634                           zdv = zh_i(ji,jl) * za_i(ji,jl)
635                           zh_i(ji,jl    ) = 0._wp
636                           za_i (ji,jl    ) = 0._wp
637                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
638                        END IF
639                     END DO
640                  ENDIF
641                  !
642               ENDIF
643               !
644               ! Compatibility tests
645               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 
646               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
647               !
648               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )
649               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
650               !
651               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ?
652               !
653               itest(4) = 1
654               DO jl = 1, i_fill
655                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0                                ! Test 4: positivity of ice concentrations
656               END DO
657               !                                         !----------------------------
658            END DO                                       ! end iteration on categories
659            !                                            !----------------------------
660         ENDIF
661      END DO
662
663      ! Add Snow in each category where za_i is not 0
664      DO jl = 1, jpl
665         DO ji = 1, idim
666            IF( za_i(ji,jl) > 0._wp ) THEN
667               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
668               ! In case snow load is in excess that would lead to transformation from snow to ice
669               ! Then, transfer the snow excess into the ice (different from icethd_dh)
670               zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 ) 
671               ! recompute h_i, h_s avoiding out of bounds values
672               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
673               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
674            ENDIF
675         END DO
676      END DO
677      !
678   END SUBROUTINE ice_var_itd
679
680
681   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
682      !!-------------------------------------------------------------------
683      !!                ***  ROUTINE ice_var_itd2   ***
684      !!
685      !! ** Purpose :  converting N-cat ice to jpl ice categories
686      !!
687      !!                  ice thickness distribution follows a gaussian law
688      !!               around the concentration of the most likely ice thickness
689      !!                           (similar as iceistate.F90)
690      !!
691      !! ** Method:   Iterative procedure
692      !!               
693      !!               1) Fill ice cat that correspond to input thicknesses
694      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
695      !!
696      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
697      !!                   by removing 25% ice area from jlmin and jlmax (resp.)
698      !!             
699      !!               3) Expand the filling to the empty cat between jlmin and jlmax
700      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
701      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
702      !!
703      !! ** Arguments : zhti: N-cat ice thickness
704      !!                zhts: N-cat snow depth
705      !!                zati: N-cat ice concentration
706      !!
707      !! ** Output    : jpl-cat
708      !!
709      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl) 
710      !!-------------------------------------------------------------------
711      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices
712      INTEGER  ::   idim, icat 
713      INTEGER, PARAMETER ::   ztrans = 0.25_wp
714      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
715      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
716      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2
717      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin
718      !!-------------------------------------------------------------------
719      !
720      idim = SIZE( zhti, 1 )
721      icat = SIZE( zhti, 2 )
722      !
723      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
724      ALLOCATE( jlmin(idim), jlmax(idim) )
725
726      ! --- initialize output fields to 0 --- !
727      zh_i(1:idim,1:jpl) = 0._wp
728      zh_s(1:idim,1:jpl) = 0._wp
729      za_i(1:idim,1:jpl) = 0._wp
730      !
731      ! --- fill the categories --- !
732      !     find where cat-input = cat-output and fill cat-output fields 
733      jlmax(:) = 0
734      jlmin(:) = 999
735      jlfil(:,:) = 0
736      DO jl1 = 1, jpl
737         DO jl2 = 1, icat
738            DO ji = 1, idim
739               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
740                  ! fill the right category
741                  zh_i(ji,jl1) = zhti(ji,jl2)
742                  zh_s(ji,jl1) = zhts(ji,jl2)
743                  za_i(ji,jl1) = zati(ji,jl2)
744                  ! record categories that are filled
745                  jlmax(ji) = MAX( jlmax(ji), jl1 )
746                  jlmin(ji) = MIN( jlmin(ji), jl1 )
747                  jlfil(ji,jl1) = jl1
748               ENDIF
749            END DO
750         END DO
751      END DO
752      !
753      ! --- fill the gaps between categories --- ! 
754      !     transfer from categories filled at the previous step to the empty ones in between
755      DO ji = 1, idim
756         jl1 = jlmin(ji)
757         jl2 = jlmax(ji)
758         IF( jl1 > 1 ) THEN
759            ! fill the lower cat (jl1-1)
760            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
761            zh_i(ji,jl1-1) = hi_mean(jl1-1)
762            ! remove from cat jl1
763            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
764         ENDIF
765         IF( jl2 < jpl ) THEN
766            ! fill the upper cat (jl2+1)
767            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
768            zh_i(ji,jl2+1) = hi_mean(jl2+1)
769            ! remove from cat jl2
770            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
771         ENDIF
772      END DO
773      !
774      jlfil2(:,:) = jlfil(:,:) 
775      ! fill categories from low to high
776      DO jl = 2, jpl-1
777         DO ji = 1, idim
778            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
779               ! fill high
780               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
781               zh_i(ji,jl) = hi_mean(jl)
782               jlfil(ji,jl) = jl
783               ! remove low
784               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
785            ENDIF
786         END DO
787      END DO
788      !
789      ! fill categories from high to low
790      DO jl = jpl-1, 2, -1
791         DO ji = 1, idim
792            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
793               ! fill low
794               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
795               zh_i(ji,jl) = hi_mean(jl) 
796               jlfil2(ji,jl) = jl
797               ! remove high
798               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
799            ENDIF
800         END DO
801      END DO
802      !
803      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
804      DEALLOCATE( jlmin, jlmax )
805      !
806   END SUBROUTINE ice_var_itd2
807
808
809   SUBROUTINE ice_var_bv
810      !!-------------------------------------------------------------------
811      !!                ***  ROUTINE ice_var_bv ***
812      !!
813      !! ** Purpose :   computes mean brine volume (%) in sea ice
814      !!
815      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
816      !!
817      !! References : Vancoppenolle et al., JGR, 2007
818      !!-------------------------------------------------------------------
819      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
820      !!-------------------------------------------------------------------
821      !
822!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
823!!   instead of setting everything to zero as just below
824      bv_i (:,:,:) = 0._wp
825      DO jl = 1, jpl
826         DO jk = 1, nlay_i
827            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
828               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
829            END WHERE
830         END DO
831      END DO
832      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
833      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
834      END WHERE
835      !
836   END SUBROUTINE ice_var_bv
837
838
839   SUBROUTINE ice_var_enthalpy
840      !!-------------------------------------------------------------------
841      !!                   ***  ROUTINE ice_var_enthalpy ***
842      !!                 
843      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
844      !!
845      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
846      !!-------------------------------------------------------------------
847      INTEGER  ::   ji, jk   ! dummy loop indices
848      REAL(wp) ::   ztmelts  ! local scalar
849      !!-------------------------------------------------------------------
850      !
851      DO jk = 1, nlay_i             ! Sea ice energy of melting
852         DO ji = 1, npti
853            ztmelts      = - tmut  * sz_i_1d(ji,jk)
854            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
855                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
856            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
857               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
858               &                    - rcp  *   ztmelts )
859         END DO
860      END DO
861      DO jk = 1, nlay_s             ! Snow energy of melting
862         DO ji = 1, npti
863            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
864         END DO
865      END DO
866      !
867   END SUBROUTINE ice_var_enthalpy
868
869#else
870   !!----------------------------------------------------------------------
871   !!   Default option         Dummy module           NO SI3 sea-ice model
872   !!----------------------------------------------------------------------
873#endif
874
875   !!======================================================================
876END MODULE icevar
Note: See TracBrowser for help on using the repository browser.