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

source: trunk/NEMO/LIM_SRC_3/limthd_lab.F90 @ 833

Last change on this file since 833 was 825, checked in by ctlod, 16 years ago

dev_002_LIM : add the LIM 3.0 component, see ticketr: #71
File size: 15.1 KB
Line 
1MODULE limthd_lab
2#if defined key_lim3
3   !!======================================================================
4   !!                       ***  MODULE limthd_lab ***
5   !!                thermodynamic lateral ablation of the ice
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   lim_thd_lab  : vertical accr./abl. and lateral ablation of sea ice
10   !! * Modules used
11   USE par_oce          ! ocean parameters
12   USE phycst           ! physical constants (OCE directory)
13   USE ice_oce          ! ice variables
14   USE thd_ice
15   USE iceini
16   USE limistate
17   USE in_out_manager
18   USE ice
19   USE par_ice
20   USE limtab
21   USE limicepoints
22     
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC lim_thd_lab        ! called by lim_thd
28
29   !! * Module variables
30   REAL(wp)  ::           &  ! constant values
31      epsi20 = 1e-20   ,  &
32      epsi13 = 1e-13   ,  &
33      zzero  = 0.e0     ,  &
34      zone   = 1.e0
35   !!----------------------------------------------------------------------
36   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2005)
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE lim_thd_lab
42       !!------------------------------------------------------------------
43       !!                ***  ROUTINE lim_thd_lab  ***
44       !!             
45       !! ** Purpose : This routine determines the time evolution of
46       !!              ice concentration, snow thickness, ice thickness
47       !!              ice temperature after lateral ablation
48       !!
49       !! ** Method  : Thermal lateral ablation
50       !!               - Heat comes from excessive bottom or total ablation
51       !!              Automatic lateral ablation
52       !!               - Concentration is "automatically reduced"
53       !!                 when it melts at the bottom
54       !!
55       !! ** Action  : 0) switches
56       !!              1) Thermal lateral ablation
57       !!                 Lead budget -> fscbq_1d, qfvbq_1d
58       !!                 Ice heat content, new ice fraction
59       !!              2) Automatic lateral ablation
60       !!              3) Constraints on lead fraction
61       !!              4) Corrections on surface and inner temperatures
62       !!              5) Variations in ice mass and volume
63       !!
64       !! References :
65       !!   Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
66       !!   Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268 
67       !!
68       !! History :
69       !!   These lines were formerly included into the vertical growth/decay routine
70       !!   original    : 01-04 (M.A. Morales Maqueda, H. Goosse, T. Fichefet)
71       !!   addition    : 02-08 (C. Ethe, G. Madec)
72       !!   addition    : 05-12 (M. Vancoppenolle ) New routine, Multilayer code,
73       !!                                           Salinity variation in the ice
74       !!------------------------------------------------------------------
75       !! * Arguments
76
77       !! * Local variables
78       INTEGER ::   ji,jj,jk,jl,index, &      ! dummy loop indices
79                    nbpb
80
81       REAL(wp) :: & 
82          zihic,      &
83          zihsn,      &
84          zihicsn,    &
85          zqice,      &
86          zqicetot,   & 
87          ziqf,       &
88          zdfrl,      &
89          zdvsnvol,   &
90          zdrfrl1,    &
91          zdfseqv
92
93       INTEGER ::     &
94          zji,        &
95          zjj
96 
97       REAL(wp), DIMENSION(jpij) ::  &
98          zfrld_old, zfrld_1d
99
100       WRITE(numout,*) 'lim_thd_lab'
101       WRITE(numout,*) '~~~~~~~~~~~'
102
103       DO jl = 1, jpl
104!--------------------------------------------------------------------------------------------------!
105! 1) Select icy points
106!--------------------------------------------------------------------------------------------------!
107         nbpb = 0
108         DO jj = 1, jpj
109            DO ji = 1, jpi
110!              IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN     
111               IF ( frld(ji,jj).lt. 1.0 ) THEN
112                  nbpb      = nbpb  + 1
113                  npb(nbpb) = (jj - 1) * jpi + ji
114               ENDIF
115            END DO
116         END DO
117
118!--------------------------------------------------------------------------------------------------!
119! 2) Convert vectors
120!--------------------------------------------------------------------------------------------------!
121
122! If there is no ice, do nothing. Otherwise, compute lateral ablation
123
124         IF ( nbpb > 0) THEN
125
126            CALL tab_2d_1d( nbpb, frld_1d    (1:nbpb)     , frld            , jpi, jpj, npb(1:nbpb) )
127            CALL tab_2d_1d( nbpb, ht_i_b     (1:nbpb)     , ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
128            CALL tab_2d_1d( nbpb, ht_s_b     (1:nbpb)     , ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
129! no use
130!           CALL tab_2d_1d( nbpb, old_ht_i_b (1:nbpb)     , old_ht_i(:,:,jl), jpi, jpj, npb(1:nbpb) )
131!           CALL tab_2d_1d( nbpb, old_ht_s_b (1:nbpb)     , old_ht_s(:,:,jl), jpi, jpj, npb(1:nbpb) )
132
133            CALL tab_2d_1d( nbpb, t_su_b     (1:nbpb)     , t_su(:,:,jl), jpi, jpj, npb(1:nbpb) )
134            CALL tab_2d_1d( nbpb, sm_i_b     (1:nbpb)     , sm_i(:,:,jl), jpi, jpj, npb(1:nbpb) )
135            DO jk = 1, nlay_s
136               CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk)     , t_s(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
137            END DO
138            DO jk = 1, nlay_i
139               CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk)     , t_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
140               CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk)     , e_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
141               CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk)     , s_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
142            END DO
143
144            CALL tab_2d_1d( nbpb, t_bo_b     (1:nbpb)     , t_bo       , jpi, jpj, npb(1:nbpb) )
145!           CALL tab_2d_1d( nbpb, thcm_1d    (1:nbpb)     , thcm       , jpi, jpj, npb(1:nbpb) )
146            CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) )
147            CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) )
148            CALL tab_2d_1d( nbpb, fseqv_1d   (1:nbpb)     , fseqv      , jpi, jpj, npb(1:nbpb) )
149
150            !dummy for temporary reasons
151            CALL tab_2d_1d( nbpb, dh_i_surf  (1:nbpb)     , dh_i_surf2D, jpi, jpj, npb(1:nbpb) )
152            CALL tab_2d_1d( nbpb, dh_i_bott  (1:nbpb)     , dh_i_bott2D, jpi, jpj, npb(1:nbpb) )
153            CALL tab_2d_1d( nbpb, fstbif_1d  (1:nbpb)     , fstbif     , jpi, jpj, npb(1:nbpb) )
154            CALL tab_2d_1d( nbpb, fsup       (1:nbpb)     , fsup2D     , jpi, jpj, npb(1:nbpb) )
155            CALL tab_2d_1d( nbpb, focea      (1:nbpb)     , focea2D    , jpi, jpj, npb(1:nbpb) )
156             
157!- Ice Sample points One-dimensional Grids (Useful only for debugging the code)
158!-----------------------------------------------------------------------------------
159      IF (ln_nicep) THEN
160
161      DO ji = 1, nbpb 
162         zji                  = MOD( npb(ji) - 1, jpi ) + 1
163         zjj                  = ( npb(ji) - 1 ) / jpi + 1
164         mask_sp_1d( ji )     = mask_sp( zji , zjj )
165         DO index = 1, arc_sp_num
166            IF ( (zji .eq.arc_sp_grid(index,1)).and.(zjj .eq.arc_sp_grid(index,2)) ) THEN
167               arc_sp_grid_1d(index) = ji*mask_sp_1d(ji) 
168            ENDIF
169            IF ( (zji .eq.ant_sp_grid(index,1)).and.(zjj .eq.ant_sp_grid(index,2)) ) THEN
170               ant_sp_grid_1d(index) = ji*mask_sp_1d(ji) 
171            ENDIF
172         END DO
173      END DO
174
175      ENDIF
176
177!-----------------------------------------------------------------------------------
178! Lateral ablation starts
179!-----------------------------------------------------------------------------------
180       DO ji = 1, nbpb
181
182          zfrld_old(ji)  = frld_1d(ji)
183          zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
184          zihsn   = 1.0 - MAX( zzero , SIGN( zone , -ht_s_b(ji) ) )
185
186          !--0) Case of total vertical ablation
187          !-------------------------------------
188          ! In the case of total ablation (all the ice ice has melted) there are only leads
189          ! frld = 1
190          frld_1d(ji)  = ( 1.0 - zihic ) + zihic * zfrld_old(ji)
191
192          !--1) Lateral ablation
193          !----------------------
194          !-- Part of solar radiation transmitted through the ice and going to the ocean
195          !-- fscbq_1d or fscmbq in 3d routines
196!         fscbq_1d(ji) = ( 1.0 - zfrld_old(ji) ) * ( 1.0 - thcm_1d(ji) ) * fstbif_1d(ji)
197          fscbq_1d(ji) = ( 1.0 - zfrld_old(ji) ) * fstbif_1d(ji)
198             !--fstbif_1d(ji) = part of the solar radiation transmitted throught the ice
199             !--thcm = portion of the solar radiation used for the lead budget...
200
201          !--Energy remaining after excessive bottom melt or total ablation
202          qfvbq_1d(ji) = fsup(ji) + ( 1.0 - zihic ) * focea(ji)
203             !!fsup latent heat remaining if bottom ablation is too strong
204             !!focea latent heat remaining in case of total ablation
205
206          !--Total lead budget
207          !!--fscbq_1d(ji) is the part of solar radiation transmitted through the ice to the ocean
208          !!               put in the lead
209          !!--qfvbq_1d is the remaining energy after excessive bottom melt or total ablation
210          !!--qldif contains -> formation of ice /snow-ice etc ...
211          !!--no heat from the ocean
212          !!--residuals from excessive bottom melt or total ablation
213          qldif_1d(ji)  = qldif_1d(ji) + qfvbq_1d(ji) + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice
214
215          !-- Total heat content per unit cell area per unit mass in the snow-ice system
216          !-- [J.kg-1.m-2]
217!         q_s_b(ji,1)     =   rcpsn*( rtt-t_s_b(ji,1) ) + xlic
218          q_s_b(ji,1)     =   rhosn*cpic*( rtt-t_s_b(ji,1) ) + rhosn*lfus
219          zqice = q_s_b(ji,1)*ht_s_b(ji) 
220          DO jk = 1, nlay_i
221             zqice = zqice + q_i_b(ji,jk)*ht_i_b(ji)/FLOAT(nlay_i)
222          END DO
223          zqicetot  = ( 1.0 - frld_1d(ji) ) * zqice
224
225          !--The concentration of ice is reduced (frld increases) if the heat
226          !--exchange between ice and ocean is positive
227         
228          ziqf = MAX( zzero , SIGN( zone ,  zqicetot - qldif_1d(ji) ) )
229          zdfrl = qldif_1d(ji) / MAX( epsi20 , zqice ) 
230          frld_1d(ji)  = ( 1.0 - ziqf )    &
231             &       +         ziqf * ( frld_1d(ji) + MAX( zzero , zdfrl ) ) 
232          !-- ice total energy per unit cell area per unit mass Q/delta_t
233          !-- fltbif_1d or ffltbif in 3d Routintes
234          fltbif_1d(ji) = ( - ( 1.0 - zfrld_old(ji) ) * zqicetot  ) / rdt_ice
235             !-->used for computing heat fluxes
236
237          !-- 2) Automatic induced lateral ablation
238          !----------------------------------------------------------------------------
239          !-- Hakkinen & Mellor, 1992.
240          zdfrl = - ( dh_i_surf(ji) + dh_i_bott(ji) ) * hakspl * ( 1.0 - zfrld_old(ji) ) &
241                / MAX( epsi13 , ht_i_b(ji) + ht_s_b(ji) * rhosn/rhoic ) 
242          zfrld_1d(ji) =  frld_1d(ji) + MAX( zzero , zdfrl )
243
244          !-- 3) Constraints on lead fraction
245          !------------------------------------------------
246          ! Sometimes, it happens that all the ice melts, but that all the heat content is not
247          ! consumed. Then, we fix the lead fraction to 0.99 (the ice concentration to 0.01)
248          ! if there is no more heat content in the ice (ziqf = 0) -> frld = 1 (only leads)
249          ! otherwise there is still some ice (ziqf = 1)           -> frld <= 0.99
250          zji = mod(npb(ji)-1,182)+1
251          zjj = ( npb(ji)-1 )/ 182 + 1
252
253          !-- 4) Update surface and internal temperature and snow/ice thicknesses
254          !-----------------------------------------------------------------------
255          t_su_b(ji)   = t_su_b(ji)   + ( 1.0 - ziqf ) * ( rtt - t_su_b(ji)   )
256          DO jk = 1, nlay_i
257             q_i_b(ji,jk) = q_i_b(ji,jk)*ziqf
258          END DO
259
260          zfrld_1d(ji) = ziqf * MIN( 0.99 * zone , zfrld_1d(ji) ) + ( 1 - ziqf )
261          zihicsn = 1.0 - MAX( zzero , SIGN( zone , - ht_i_b(ji) - ht_s_b(ji) ) )
262          zfrld_1d(ji) = ( 1.0 - zihicsn ) + zihicsn * zfrld_1d(ji)
263   
264          !-- 5) Variation of ice/snow volume and mass
265          !-------------------------------------------
266          dvlbq_1d(ji)   = zihic * ( zfrld_old(ji) - frld_1d(ji) ) * ht_i_b(ji)
267          rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic
268          zdvsnvol    = zihsn * ( zfrld_old(ji) - frld_1d(ji) ) * ht_s_b(ji)
269          rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn
270
271!        Equivalent salt flux (3) Lateral Ice melt
272!        -----------------------------------------
273         fseqv_1d(ji)   = fseqv_1d(ji)   + &
274                          ( soce - sm_i_b(ji) ) * &
275                          MIN( dvlbq_1d(ji) , 0.0) * rhoic / rdt_ice
276
277          !-- 6) Final update
278          !------------------
279          ht_s_b(ji)  = ziqf * ht_s_b(ji)
280          zdrfrl1 = ziqf * ( 1.0 -  frld_1d(ji) ) / MAX( epsi20 , 1.0 - zfrld_1d(ji) )
281          ht_s_b(ji) = zdrfrl1 * ht_s_b(ji)
282          ht_i_b(ji) = zdrfrl1 * ht_i_b(ji)
283          frld_1d(ji)    = zfrld_1d(ji)
284
285       END DO
286
287!-4.4) Back to the geographic grid
288!----------------------------------------------
289         CALL tab_1d_2d( nbpb, frld       , npb, frld_1d(1:nbpb), jpi, jpj )
290         CALL tab_1d_2d( nbpb, ht_i(:,:,jl), npb, ht_i_b(1:nbpb), jpi, jpj )
291         CALL tab_1d_2d( nbpb, ht_s(:,:,jl), npb, ht_s_b(1:nbpb), jpi, jpj )
292         CALL tab_1d_2d( nbpb, t_su(:,:,jl), npb, t_su_b(1:nbpb), jpi, jpj )
293         CALL tab_1d_2d( nbpb, sm_i(:,:,jl), npb, sm_i_b(1:nbpb), jpi, jpj )
294         DO jk = 1, nlay_s
295            CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b(1:nbpb,jk), jpi, jpj)
296         END DO
297         DO jk = 1, nlay_i
298            CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b(1:nbpb,jk), jpi, jpj)
299            CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b(1:nbpb,jk), jpi, jpj)
300            CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b(1:nbpb,jk), jpi, jpj)
301         END DO
302
303         CALL tab_1d_2d( nbpb, qldif      , npb, qldif_1d  (1:nbpb)     , jpi, jpj )
304         CALL tab_1d_2d( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj )
305
306         CALL tab_1d_2d( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj )
307         CALL tab_1d_2d( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj )
308         CALL tab_1d_2d( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj )
309
310      ENDIF
311       
312!-4.5) Update sea ice thickness   
313!----------------------------------------------------------
314      DO jj = 1, jpj
315         DO ji = 1, jpi
316            phicif(ji,jj) = ht_i(ji,jj,1) 
317!           ht_i(ji,jj,1) = ht_i(ji,jj,1) *  ( zone -  MAX( zzero, SIGN( zone, - at_i(ji,jj) ) ) )
318            ht_i(ji,jj,1) = ht_i(ji,jj,1) *  ( zone -  MAX( zzero, SIGN( zone, - ( 1.0 - frld(ji,jj) ) ) ) ) 
319         END DO
320      END DO
321
322     
323    END DO
324    END SUBROUTINE lim_thd_lab
325#else
326   !!======================================================================
327   !!                       ***  MODULE limthd_lab   ***
328   !!                              no sea ice model 
329   !!======================================================================
330CONTAINS
331   SUBROUTINE lim_thd_lab          ! Empty routine
332   END SUBROUTINE lim_thd_lab
333#endif
334 END MODULE limthd_lab
Note: See TracBrowser for help on using the repository browser.