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.
limcat_1D.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90 @ 4765

Last change on this file since 4765 was 4333, checked in by clem, 10 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config

File size: 9.5 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, zdh   ! 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         
88         IF( zhti(ji) > 0._wp ) THEN
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               ! *** case ice is thicker: fill categories >1
111            ELSE
112
113               ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2
114               DO jl = 1, i_fill - 1
115                  zht_i(ji,jl) = ( zhi_max(jl) + zhi_max(jl-1) ) * 0.5_wp
116               END DO
117               
118               ! find which category (jl0) the input ice thickness falls into
119               jl0 = i_fill
120               DO jl = 1, i_fill
121                  IF ( ( zhti(ji) >= zhi_max(jl-1) ) .AND. ( zhti(ji) < zhi_max(jl) ) ) THEN
122                     jl0 = jl
123           CYCLE
124                  ENDIF
125               END DO
126               
127               ! Concentrations in the (i_fill-1) categories
128               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
129               DO jl = 1, i_fill - 1
130                  IF ( jl == jl0 ) CYCLE
131                  zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
132                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
133               END DO
134               
135               ! Concentration in the last (i_fill) category
136               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
137               
138               ! Ice thickness in the last (i_fill) category
139               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
140               zht_i(ji,i_fill) = ( zhti(ji)*zai(ji) -  zV ) / za_i(ji,i_fill) 
141               
142            ENDIF ! case ice is thick or thin
143           
144            !---------------------
145            ! Compatibility tests
146            !---------------------
147            ! Test 1: area conservation
148            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
149            IF ( zconv < epsi06 ) ztest_1 = 1
150           
151            ! Test 2: volume conservation
152            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
153            IF ( zconv < epsi06 ) ztest_2 = 1
154           
155            ! Test 3: thickness of the last category is in-bounds ?
156            IF ( zht_i(ji,i_fill) >= zhi_max(i_fill-1) ) ztest_3 = 1
157           
158            ! Test 4: positivity of ice concentrations
159            ztest_4 = 1
160            DO jl = 1, i_fill
161               IF ( za_i(ji,jl) < 0._wp ) ztest_4 = 0
162            END DO
163           
164            ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
165                                                           !============================
166         END DO                                            ! end iteration on categories
167                                                           !============================
168         ! Check if tests have passed (i.e. volume conservation...)
169         !IF ( ztests /= 4 ) THEN
170         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
171         !   WRITE(numout,*) ' !! ALERT categories distribution !!'
172         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
173         !   WRITE(numout,*) ' *** ztests is not equal to 4 '
174         !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
175         !   WRITE(numout,*) 'i_fill=',i_fill
176         !   WRITE(numout,*) 'zai(ji)=',zai(ji)
177         !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)
178         !ENDIF
179         
180         ENDIF ! if zhti > 0
181
182      END DO ! i loop
183
184      ! ------------------------------------------------
185      ! Adding Snow in each category where za_i is not 0
186      ! ------------------------------------------------
187      DO jl = 1, jpl
188         DO ji = 1, ijpij
189            IF( za_i(ji,jl) > 0._wp ) THEN
190               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
191               ! In case snow load is in excess that would lead to transformation from snow to ice
192               ! Then, transfer the snow excess into the ice (different from limthd_dh)
193               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
194               ! recompute ht_i, ht_s avoiding out of bounds values
195               zht_i(ji,jl) = MIN( zhi_max(jl), zht_i(ji,jl) + zdh )
196               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
197            ENDIF
198         ENDDO
199      ENDDO
200
201      IF( nn_timing == 1 )  CALL timing_stop('limcat_1D')
202     
203    END SUBROUTINE lim_cat_1D
204   
205#else
206   !!----------------------------------------------------------------------
207   !!   Default option :         Empty module          NO LIM sea-ice model
208   !!----------------------------------------------------------------------
209CONTAINS
210   SUBROUTINE lim_cat_1D          ! Empty routine
211   END SUBROUTINE lim_cat_1D
212#endif
213
214   !!======================================================================
215END MODULE limcat_1D
Note: See TracBrowser for help on using the repository browser.