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_lac_2.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released

  • Property svn:keywords set to Id
File size: 12.7 KB
Line 
1MODULE limthd_lac_2
2#if defined key_lim2
3   !!======================================================================
4   !!                       ***  MODULE limthd_lac_2   ***
5   !!                lateral thermodynamic growth of the ice
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   lim_lat_acr_2    : lateral accretion of ice
10   !! * Modules used
11   USE par_oce          ! ocean parameters
12   USE phycst
13   USE thd_ice_2
14   USE ice_2
15   USE limistate_2 
16     
17   IMPLICIT NONE
18   PRIVATE
19
20   !! * Routine accessibility
21   PUBLIC lim_thd_lac_2   ! called by lim_thd_2
22
23   !! * Module variables
24   REAL(wp)  ::           &  ! constant values
25      epsi20 = 1.e-20  ,  &
26      epsi13 = 1.e-13  ,  &
27      zzero  = 0.e0    ,  &
28      zone   = 1.e0
29   !!----------------------------------------------------------------------
30   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
31   !! $Id$
32   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34CONTAINS
35   
36   SUBROUTINE lim_thd_lac_2( kideb, kiut )
37      !!-------------------------------------------------------------------
38      !!               ***   ROUTINE lim_thd_lac_2  ***
39      !! 
40      !! ** Purpose : Computation of the evolution of the ice thickness and
41      !!      concentration as a function of the heat balance in the leads.
42      !!      It is only used for lateral accretion
43      !!       
44      !! ** Method  : Ice is formed in the open water when ocean lose heat
45      !!      (heat budget of open water Bl is negative) .
46      !!      Computation of the increase of 1-A (ice concentration) fol-
47      !!      lowing the law :
48      !!      (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
49      !!       where - h0 is the thickness of ice created in the lead
50      !!             - a is a minimum fraction for leads
51      !!             - F is a monotonic non-increasing function defined as:
52      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
53      !!             - exld is the exponent closure rate (=2 default val.)
54      !!
55      !! ** Action : - Adjustment of snow and ice thicknesses and heat
56      !!                content in brine pockets
57      !!             - Updating ice internal temperature
58      !!             - Computation of variation of ice volume and mass
59      !!             - Computation of frldb after lateral accretion and
60      !!               update h_snow_1d, h_ice_1d and tbif_1d(:,:)     
61      !!
62      !! ** References :
63      !!      M. Maqueda, 1995, PhD Thesis, Univesidad Complutense Madrid
64      !!      Fichefet T. and M. Maqueda 1997, J. Geo. Res., 102(C6),
65      !!                                                12609 -12646   
66      !! History :
67      !!   1.0  !  01-04 (LIM)  original code
68      !!   2.0  !  02-08 (C. Ethe, G. Madec)  F90, mpp
69      !!-------------------------------------------------------------------
70      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
71      USE wrk_nemo, ONLY: wrk_1d_1, wrk_1d_2, wrk_1d_3, wrk_1d_4, wrk_1d_5, wrk_1d_6
72      USE in_out_manager, ONLY: ctl_stop
73      !!
74      !! * Arguments
75      INTEGER , INTENT(IN)::  &
76         kideb          ,   &  ! start point on which the the computation is applied
77         kiut                  ! end point on which the the computation is applied
78
79      !! * Local variables
80      INTEGER ::            &
81         ji             ,   &  !  dummy loop indices
82         iicefr         ,   &  !  1 = existing ice ; 0 = no ice
83         iiceform       ,   &  !  1 = ice formed   ; 0 = no ice formed
84         ihemis                !  dummy indice
85      REAL(wp), POINTER, DIMENSION(:) :: &
86         zqbgow           ,  &  !  heat budget of the open water (negative)
87         zfrl_old         ,  &  !  previous sea/ice fraction
88         zhice_old        ,  &  !  previous ice thickness
89         zhice0           ,  &  !  thickness of newly formed ice in leads
90         zfrlmin          ,  &  !  minimum fraction for leads
91         zdhicbot               !  part of thickness of newly formed ice in leads which
92                                !  has been already used in transport for example
93      REAL(wp)  ::  &
94         zhemis           ,  &  !  hemisphere (0 = North, 1 = South)
95         zhicenew         ,  &  !  new ice thickness
96         zholds2          ,  &  !  ratio of previous ice thickness and 2
97         zhnews2          ,  &  !  ratio of new ice thickness and 2
98         zfrlnew          ,  &  !  new sea/ice fraction
99         zfrld            ,  &  !  ratio of sea/ice fraction and minimum fraction for leads
100         zfrrate          ,  &  !  leads-closure rate
101         zdfrl                  !  sea-ice fraction increment
102      REAL(wp)  ::  &
103         zdh1 , zdh2 , zdh3 , zdh4, zdh5   , &   ! tempory scalars
104         ztint , zta1 , zta2 , zta3 , zta4 , &
105         zah, zalpha , zbeta
106      !!---------------------------------------------------------------------     
107               
108      IF(wrk_in_use(1, 1,2,3,4,5,6))THEN
109         CALL ctl_stop('lim_thd_lac_2 : requestead workspace arrays unavailable.')
110         RETURN
111      END IF
112      ! Set-up pointers to sub-arrays of workspace arrays
113      zqbgow    => wrk_1d_1(1:jpij)
114      zfrl_old  => wrk_1d_2(1:jpij)         
115      zhice_old => wrk_1d_3(1:jpij)       
116      zhice0    => wrk_1d_4(1:jpij)       
117      zfrlmin   => wrk_1d_5(1:jpij)       
118      zdhicbot  => wrk_1d_6(1:jpij) 
119     
120      !--------------------------------------------------------------
121      !   Computation of the heat budget of the open water (negative)
122      !--------------------------------------------------------------
123     
124      DO ji = kideb , kiut     
125         zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)
126      END DO
127     
128      !-----------------------------------------------------------------
129      !   Taking the appropriate values for the corresponding hemisphere
130      !-----------------------------------------------------------------
131      DO ji = kideb , kiut
132         zhemis       = MAX( zzero , SIGN( zone , frld_1d(ji) - 2.0 ) ) 
133         ihemis       = INT( 1 + zhemis )
134         zhice0  (ji) = hiccrit( ihemis ) 
135         zfrlmin (ji) = acrit  ( ihemis )   
136         frld_1d (ji) = frld_1d(ji) - 2.0 * zhemis
137         zfrl_old(ji) = frld_1d(ji)
138      END DO
139     
140      !-------------------------------------------------------------------
141      !     Lateral Accretion (modification of the fraction of open water)
142      !     The ice formed in the leads has always a thickness zhice0, but
143      !     only a fraction zfrrate of the ice formed contributes to the
144      !     increase of the ice fraction. The remaining part (1-zfrrate)
145      !     is rather assumed to lead to an increase in the thickness of the
146      !     pre-existing ice (transport for example).
147      !     Morales Maqueda, 1995 - Fichefet and Morales Maqueda, 1997
148      !---------------------------------------------------------------------
149     
150!CDIR NOVERRCHK
151      DO ji = kideb , kiut
152         iicefr       = 1 - MAX( 0, INT( SIGN( 1.5 * zone , zfrl_old(ji) - 1.0 + epsi13 ) ) )
153         !---computation of the leads-closure rate
154         zfrld        = MIN( zone , ( 1.0 - frld_1d(ji) ) / ( 1.0 - zfrlmin(ji) ) )
155         zfrrate      = ( 1.0 - zfrld**exld )**( 1.0 / exld )
156         !--computation of the sea-ice fraction increment and the new fraction
157         zdfrl        = ( zfrrate / zhice0(ji) )  * ( zqbgow(ji) / xlic )
158         zfrlnew      = zfrl_old(ji) + zdfrl
159         !--update the sea-ice fraction
160         frld_1d   (ji) = MAX( zfrlnew , zfrlmin(ji) )
161         !--computation of the remaining part of ice thickness which has been already used
162         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) ) & 
163                      -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic ) 
164      END DO
165 
166      !----------------------------------------------------------------------------------------
167      !      Ajustement of snow and ice thicknesses and updating the total heat stored in brine pockets 
168      !      The thickness of newly formed ice is averaged with that of the pre-existing
169      !         (1-Anew) * hinew = (1-Aold) * hiold + ((1-Anew)-(1-Aold)) * h0
170      !      Snow is distributed over the new ice-covered area
171      !         (1-Anew) * hsnew = (1-Aold) * hsold           
172      !--------------------------------------------------------------------------------------------
173     
174      DO ji = kideb , kiut
175         iicefr       = 1 - MAX( 0, INT( SIGN( 1.5 * zone , zfrl_old(ji) - 1.0 + epsi13 ) ) )
176         zhice_old(ji) = h_ice_1d(ji)
177         zhicenew      = iicefr * zhice_old(ji) + ( 1 - iicefr ) * zhice0(ji)
178         zalpha        = ( 1. - zfrl_old(ji) ) / ( 1.- frld_1d(ji) )
179         h_snow_1d(ji) = zalpha * h_snow_1d(ji)
180         h_ice_1d (ji) = zalpha * zhicenew + ( 1.0 - zalpha ) * zhice0(ji)
181         qstbif_1d(ji) = zalpha * qstbif_1d(ji) 
182      END DO
183     
184      !-------------------------------------------------------
185      !   Ajustement of ice internal temperatures
186      !-------------------------------------------------------
187     
188      DO ji = kideb , kiut
189         iicefr      = 1 - MAX( 0, INT( SIGN( 1.5 * zone , zfrl_old(ji) - 1.0 + epsi13 ) ) )
190         iiceform    = 1 - MAX( 0 ,INT( SIGN( 1.5 * zone , zhice0(ji) - h_ice_1d(ji) ) ) )
191         zholds2     = zhice_old(ji)/ 2.
192         zhnews2     = h_ice_1d(ji) / 2.
193         zdh1        = MAX( zzero ,  zhice_old(ji)   - zhnews2 )
194         zdh2        = MAX( zzero , -zhice_old(ji)   + zhnews2 )
195         zdh3        = MAX( zzero ,  h_ice_1d(ji) - zholds2 )
196         zdh4        = MAX( zzero , -h_ice_1d(ji) + zholds2 )
197         zdh5        = MAX( zzero , zhice0(ji)      - zholds2 )
198         ztint       =       iiceform   * (  ( zholds2 - zdh3 ) * tbif_1d(ji,3) + zdh4 * tbif_1d(ji,2) )      &
199            &                           / MAX( epsi20 , h_ice_1d(ji) - zhice0(ji) )                           &
200            &                 + ( 1 - iiceform ) * tfu_1d(ji)
201         zta1        = iicefr * ( 1.  - zfrl_old(ji) ) * tbif_1d(ji,2) 
202         zta2        = iicefr * ( 1.  - zfrl_old(ji) ) * tbif_1d(ji,3)
203         zta3        = iicefr * ( 1.  - zfrl_old(ji) ) * ztint
204         zta4        = ( zfrl_old(ji) - frld_1d   (ji) ) * tfu_1d(ji)
205         zah         = ( 1. - frld_1d(ji) ) * zhnews2 
206
207         tbif_1d(ji,2) = (  MIN( zhnews2 , zholds2 )                                              * zta1   &
208            &          + ( 1 - iiceform ) * ( zholds2 - zdh1 )                                    * zta2   &
209            &          + ( iiceform * ( zhnews2 - zhice0(ji) + zdh5 ) + ( 1 - iiceform ) * zdh2 ) * zta3   & 
210            &          + MIN ( zhnews2 , zhice0(ji) )                                             * zta4   &
211            &          ) / zah
212         
213         tbif_1d(ji,3) =     (  iiceform * ( zhnews2 - zdh3 )                                          * zta1  &
214            &              + ( iiceform * zdh3 + ( 1 - iiceform ) * zdh1 )                             * zta2  &
215            &              + ( iiceform * ( zhnews2 - zdh5 ) + ( 1 - iiceform ) * ( zhnews2 - zdh1 ) ) * zta3  & 
216            &              + ( iiceform * zdh5 + ( 1 - iiceform ) * zhnews2 )                          * zta4  &
217            &            ) / zah
218         !---removing the remaining part of ice formed which has been already used
219         zbeta         = h_ice_1d(ji) / ( h_ice_1d(ji) + zdhicbot(ji) )
220         h_ice_1d(ji)  = h_ice_1d(ji) + zdhicbot(ji)
221         tbif_1d (ji,2)= zbeta * tbif_1d(ji,2) + ( 1.0 - zbeta ) * tbif_1d(ji,3)
222         tbif_1d (ji,3)= ( 2. * zbeta - 1.0 ) * tbif_1d(ji,3) + ( 2. * zdhicbot(ji) / h_ice_1d(ji) ) * tfu_1d(ji)
223         
224      END DO
225     
226      !-------------------------------------------------------------
227      !    Computation of variation of ice volume and ice mass
228      !           Vold = (1-Aold) * hiold ; Vnew = (1-Anew) * hinew
229      !           dV = Vnew - Vold
230      !-------------------------------------------------------------
231     
232      DO ji = kideb , kiut
233         dvlbq_1d  (ji) = ( 1. - frld_1d(ji) ) * h_ice_1d(ji) - ( 1. - zfrl_old(ji) ) * zhice_old(ji)
234         rdmicif_1d(ji) = rdmicif_1d(ji) + rhoic * dvlbq_1d(ji)
235      END DO
236     
237      IF(wrk_not_released(1, 1,2,3,4,5,6))THEN
238         CALL ctl_stop('lim_thd_lac_2 : failed to release workspace arrays.')
239      END IF
240
241   END SUBROUTINE lim_thd_lac_2
242#else
243   !!======================================================================
244   !!                       ***  MODULE limthd_lac_2   ***
245   !!                           no sea ice model
246   !!======================================================================
247CONTAINS
248   SUBROUTINE lim_thd_lac_2           ! Empty routine
249   END SUBROUTINE lim_thd_lac_2
250#endif
251END MODULE limthd_lac_2
Note: See TracBrowser for help on using the repository browser.