source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90 @ 4298

Last change on this file since 4298 was 4298, checked in by clem, 7 years ago
File size: 8.7 KB
Line 
1MODULE limcat_1D
2   !!======================================================================
3   !!                     ***  MODULE  limcat_1D  ***
4   !!  Used for LIM3 to convert cell averages of ice thickness, snow thickness
5   !!  and ice cover into a prescribed distribution over the cell.
6   !!  (Example of application: BDY forcings when input are cell averaged) 
7   !!======================================================================
8   !! History :   -   ! Original code from M. Vancoppenolle (?)
9   !!                 ! 2011-12 (C. Rousset) rewritten for clarity
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3' :                                    LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_cat_1D      :  main subroutine
16   !!----------------------------------------------------------------------
17   !! Modules used
18   USE phycst
19   USE oce             ! dynamics and tracers variables
20   USE dom_oce
21   USE sbc_oce         ! Surface boundary condition: ocean fields
22   USE par_ice         ! ice parameters
23   USE ice             ! ice variables
24   USE eosbn2          ! equation of state
25   USE in_out_manager
26   USE dom_ice
27   USE ice
28   USE lbclnk
29   USE timing          ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! Accessibility
35   PUBLIC lim_cat_1D     
36
37CONTAINS
38
39   SUBROUTINE lim_cat_1D(zhti,zhts,zai,zht_i,zht_s,za_i)
40      !! Local variables
41      INTEGER  :: ji, jk, jl             ! dummy loop indices
42      INTEGER  :: ijpij, i_fill, jl0, ztest_1, ztest_2, ztest_3, ztest_4, ztests 
43      REAL(wp) :: zarg, zV, zconv
44      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
45      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
46      REAL(wp) :: epsi06 = 1.0e-6
47      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
48      REAL(wp), DIMENSION(0:jpl)   ::   zhi_max         !:Boundary of ice thickness categories in thickness space
49 
50      IF( nn_timing == 1 )  CALL timing_start('limcat_1D')
51      !--------------------------------------------------------------------
52      ! initialisation of variables
53      !--------------------------------------------------------------------
54      ijpij = SIZE(zhti,1)
55      zht_i(1:ijpij,1:jpl) = 0._wp
56      zht_s(1:ijpij,1:jpl) = 0._wp
57      za_i (1:ijpij,1:jpl) = 0._wp
58
59      !------------------------------------------------------------------------------------
60      ! Distribute ice concentration and thickness into the categories
61      !------------------------------------------------------------------------------------
62      ! Method: we first try to fill the jpl ice categories bounded by thicknesses
63      !         hmax(0:jpl) with a gaussian distribution, and check whether the distribution
64      !         fulfills volume and area conservation, positivity and ice categories bounds.
65      !         In other words, if ice input is too thin, the last category (jpl)
66      !         cannot be filled, so we try to fill jpl-1 categories...
67      !         And so forth iteratively until the number of categories filled
68      !         fulfills ice volume concervation between input and output (ztests=4)
69      !--------------------------------------------------------------------------------------
70
71      !- Thickness categories boundaries
72      !    hi_max is calculated in iceini.F90 but since limcat_1D.F90 routine
73      !    is called before (in bdydta.F90), one must recalculate it.   
74      !    Note clem: there may be a way of doing things cleaner
75      !----------------------------------
76      zhi_max(:) = 0._wp
77      zc1 =  3._wp / REAL( jpl , wp ) ; zc2 = 10._wp * zc1 ; zc3 =  3._wp
78      DO jl = 1, jpl
79         zx1 = REAL( jl-1 , wp ) / REAL( jpl , wp )
80         zhi_max(jl) = zhi_max(jl-1) + zc1 + zc2 * ( 1._wp + TANH( zc3 * ( zx1 - 1._wp ) ) )
81      END DO
82 
83      ! ----------------------------------------
84      ! distribution over the jpl ice categories
85      ! ----------------------------------------
86      DO ji = 1, ijpij
87        ! snow thickness in each category
88         zht_s(ji,1:jpl) = zhts(ji)
89         
90         ! initialisation of tests
91         ztest_1 = 0
92         ztest_2 = 0
93         ztest_3 = 0
94         ztest_4 = 0
95         ztests  = 0
96         
97         i_fill = jpl + 1                                    !====================================
98         DO WHILE ( ( ztests /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
99            ! iteration                                      !====================================
100            i_fill = i_fill - 1
101           
102            ! initialisation of ice variables for each try
103            zht_i(ji,1:jpl) = 0._wp
104            za_i (ji,1:jpl) = 0._wp
105           
106            ! *** case very thin ice: fill only category 1
107            IF ( i_fill == 1 ) THEN
108               zht_i(ji,1) = zhti(ji)
109               za_i (ji,1) = zai (ji)
110               
111               ! *** case ice is thicker: fill categories >1
112            ELSE
113
114               ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2
115               DO jl = 1, i_fill - 1
116                  zht_i(ji,jl) = ( zhi_max(jl) + zhi_max(jl-1) ) * 0.5_wp
117               END DO
118               
119               ! find which category (jl0) the input ice thickness falls into
120               jl0 = i_fill
121               DO jl = 1, i_fill
122                  IF ( ( zhti(ji) >= zhi_max(jl-1) ) .AND. ( zhti(ji) < zhi_max(jl) ) ) THEN
123                     jl0 = jl
124           CYCLE
125                  ENDIF
126               END DO
127               
128               ! Concentrations in the (i_fill-1) categories
129               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
130               DO jl = 1, i_fill - 1
131                  IF ( jl == jl0 ) CYCLE
132                  zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
133                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
134               END DO
135               
136               ! Concentration in the last (i_fill) category
137               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
138               
139               ! Ice thickness in the last (i_fill) category
140               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
141               zht_i(ji,i_fill) = ( zhti(ji)*zai(ji) -  zV ) / za_i(ji,i_fill) 
142               
143            ENDIF ! case ice is thick or thin
144           
145            !---------------------
146            ! Compatibility tests
147            !---------------------
148            ! Test 1: area conservation
149            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
150            IF ( zconv < epsi06 ) ztest_1 = 1
151           
152            ! Test 2: volume conservation
153            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
154            IF ( zconv < epsi06 ) ztest_2 = 1
155           
156            ! Test 3: thickness of the last category is in-bounds ?
157            IF ( zht_i(ji,i_fill) >= zhi_max(i_fill-1) ) ztest_3 = 1
158           
159            ! Test 4: positivity of ice concentrations
160            ztest_4 = 1
161            DO jl = 1, i_fill
162               IF ( za_i(ji,jl) < 0._wp ) ztest_4 = 0
163            END DO
164           
165            ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
166                                                           !============================
167         END DO                                            ! end iteration on categories
168                                                           !============================
169         ! Check if tests have passed (i.e. volume conservation...)
170         !IF ( ztests /= 4 ) THEN
171         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
172         !   WRITE(numout,*) ' !! ALERT categories distribution !!'
173         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
174         !   WRITE(numout,*) ' *** ztests is not equal to 4 '
175         !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
176         !   WRITE(numout,*) 'i_fill=',i_fill
177         !   WRITE(numout,*) 'zai(ji)=',zai(ji)
178         !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)
179         !ENDIF
180           
181      END DO ! i loop
182
183      IF( nn_timing == 1 )  CALL timing_stop('limcat_1D')
184     
185    END SUBROUTINE lim_cat_1D
186   
187#else
188   !!----------------------------------------------------------------------
189   !!   Default option :         Empty module          NO LIM sea-ice model
190   !!----------------------------------------------------------------------
191CONTAINS
192   SUBROUTINE lim_cat_1D          ! Empty routine
193   END SUBROUTINE lim_cat_1D
194#endif
195
196   !!======================================================================
197END MODULE limcat_1D
Note: See TracBrowser for help on using the repository browser.