source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90 @ 4161

Last change on this file since 4161 was 4161, checked in by cetlod, 8 years ago

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

File size: 8.0 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 
48      IF( nn_timing == 1 )  CALL timing_start('limcat_1D')
49      !--------------------------------------------------------------------
50      ! initialisation of variables
51      !--------------------------------------------------------------------
52      ijpij = SIZE(zhti,1)
53      zht_i(1:ijpij,1:jpl) = 0.d0
54      zht_s(1:ijpij,1:jpl) = 0.d0
55      za_i (1:ijpij,1:jpl) = 0.d0
56
57      !------------------------------------------------------------------------------------
58      ! Distribute ice concentration and thickness into the categories
59      !------------------------------------------------------------------------------------
60      ! Method: we first try to fill the jpl ice categories bounded by thicknesses
61      !         hmax(0:jpl) with a gaussian distribution, and check whether the distribution
62      !         fulfills volume and area conservation, positivity and ice categories bounds.
63      !         In other words, if ice input is too thin, the last category (jpl)
64      !         cannot be filled, so we try to fill jpl-1 categories...
65      !         And so forth iteratively until the number of categories filled
66      !         fulfills ice volume concervation between input and output (ztests=4)
67      !--------------------------------------------------------------------------------------
68 
69      ! ----------------------------------------
70      ! distribution over the jpl ice categories
71      ! ----------------------------------------
72      DO ji = 1, ijpij
73        ! snow thickness in each category
74         zht_s(ji,1:jpl) = zhts(ji)
75         
76         ! initialisation of tests
77         ztest_1 = 0
78         ztest_2 = 0
79         ztest_3 = 0
80         ztest_4 = 0
81         ztests  = 0
82         
83         i_fill = jpl + 1                                    !====================================
84         DO WHILE ( ( ztests /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories 
85            ! iteration                                      !====================================
86            i_fill = i_fill - 1
87           
88            ! initialisation of ice variables for each try
89            zht_i(ji,1:jpl) = 0.d0
90            za_i (ji,1:jpl) = 0.d0
91           
92            ! *** case very thin ice: fill only category 1
93            IF ( i_fill == 1 ) THEN
94               zht_i(ji,1) = zhti(ji)
95               za_i (ji,1) = zai (ji)
96               
97               ! *** case ice is thicker: fill categories >1
98            ELSE
99
100               ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2
101               DO jl = 1, i_fill - 1
102                  zht_i(ji,jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2.
103               END DO
104               
105               ! find which category (jl0) the input ice thickness falls into
106               jl0 = i_fill
107               DO jl = 1, i_fill
108                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
109                     jl0 = jl
110           CYCLE
111                  ENDIF
112               END DO
113               
114               ! Concentrations in the (i_fill-1) categories
115               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
116               DO jl = 1, i_fill - 1
117                  IF ( jl == jl0 ) CYCLE
118                  zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) / 2. )
119                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
120               END DO
121               
122               ! Concentration in the last (i_fill) category
123               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
124               
125               ! Ice thickness in the last (i_fill) category
126               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
127               zht_i(ji,i_fill) = ( zhti(ji)*zai(ji) -  zV ) / za_i(ji,i_fill) 
128               
129            ENDIF ! case ice is thick or thin
130           
131            !---------------------
132            ! Compatibility tests
133            !---------------------
134            ! Test 1: area conservation
135            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
136            IF ( zconv < epsi06 ) ztest_1 = 1
137           
138            ! Test 2: volume conservation
139            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
140            IF ( zconv < epsi06 ) ztest_2 = 1
141           
142            ! Test 3: thickness of the last category is in-bounds ?
143            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) ztest_3 = 1
144           
145            ! Test 4: positivity of ice concentrations
146            ztest_4 = 1
147            DO jl = 1, i_fill
148               IF ( za_i(ji,jl) < 0.0d0 ) ztest_4 = 0
149            END DO
150           
151            ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
152                                                           !============================
153         END DO                                            ! end iteration on categories
154                                                           !============================
155         ! Check if tests have passed (i.e. volume conservation...)
156         !IF ( ztests /= 4 ) THEN
157         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
158         !   WRITE(numout,*) ' !! ALERT categories distribution !!'
159         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
160         !   WRITE(numout,*) ' *** ztests is not equal to 4 '
161         !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
162         !   WRITE(numout,*) 'i_fill=',i_fill
163         !   WRITE(numout,*) 'zai(ji)=',zai(ji)
164         !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)
165         !ENDIF
166           
167      END DO ! i loop
168
169      IF( nn_timing == 1 )  CALL timing_stop('limcat_1D')
170     
171    END SUBROUTINE lim_cat_1D
172   
173#else
174   !!----------------------------------------------------------------------
175   !!   Default option :         Empty module          NO LIM sea-ice model
176   !!----------------------------------------------------------------------
177CONTAINS
178   SUBROUTINE lim_cat_1D          ! Empty routine
179   END SUBROUTINE lim_cat_1D
180#endif
181
182   !!======================================================================
183END MODULE limcat_1D
Note: See TracBrowser for help on using the repository browser.