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 branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90 @ 4345

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

clean the glob_sum diags and make scalar.nc work

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.