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.
limistate.F90 in branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 6347

Last change on this file since 6347 was 6347, checked in by gm, 8 years ago

#1683: SIMPLIF-1 : Phase with the v3.6_Stable (DOC+ZDF+traqsr+lbedo)

  • Property svn:keywords set to Id
File size: 28.8 KB
Line 
1MODULE limistate
2   !!======================================================================
3   !!                     ***  MODULE  limistate  ***
4   !!              Initialisation of diagnostics ice variables
5   !!======================================================================
6   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
7   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
8   !!             -   ! 2014    (C. Rousset) add N/S initializations
9   !!----------------------------------------------------------------------
10#if defined key_lim3
11   !!----------------------------------------------------------------------
12   !!   'key_lim3' :                                    LIM3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   lim_istate      :  Initialisation of diagnostics ice variables
15   !!   lim_istate_init :  initialization of ice state and namelist read
16   !!----------------------------------------------------------------------
17   USE phycst           ! physical constant
18   USE oce              ! dynamics and tracers variables
19   USE dom_oce          ! ocean domain
20   USE sbc_oce          ! Surface boundary condition: ocean fields
21   USE sbc_ice          ! Surface boundary condition: ice fields
22   USE eosbn2           ! equation of state
23   USE ice              ! sea-ice variables
24   USE par_oce          ! ocean parameters
25   USE dom_ice          ! sea-ice domain
26   USE in_out_manager   ! I/O manager
27   USE lib_mpp          ! MPP library
28   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
29   USE wrk_nemo         ! work arrays
30   USE fldread          ! read input fields
31   USE iom
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   lim_istate      ! routine called by lim_init.F90
37
38   !                          !!** init namelist (namiceini) **
39   REAL(wp) ::   rn_thres_sst   ! threshold water temperature for initial sea ice
40   REAL(wp) ::   rn_hts_ini_n   ! initial snow thickness in the north
41   REAL(wp) ::   rn_hts_ini_s   ! initial snow thickness in the south
42   REAL(wp) ::   rn_hti_ini_n   ! initial ice thickness in the north
43   REAL(wp) ::   rn_hti_ini_s   ! initial ice thickness in the south
44   REAL(wp) ::   rn_ati_ini_n   ! initial leads area in the north
45   REAL(wp) ::   rn_ati_ini_s   ! initial leads area in the south
46   REAL(wp) ::   rn_smi_ini_n   ! initial salinity
47   REAL(wp) ::   rn_smi_ini_s   ! initial salinity
48   REAL(wp) ::   rn_tmi_ini_n   ! initial temperature
49   REAL(wp) ::   rn_tmi_ini_s   ! initial temperature
50
51   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
52   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point
53   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point
54   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
55   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
56   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
57   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si    ! structure of input fields (file informations, fields read)
59
60   LOGICAL  ::  ln_iceini        ! initialization or not
61   LOGICAL  ::  ln_iceini_file   ! Ice initialization state from 2D netcdf file
62   !!----------------------------------------------------------------------
63   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
64   !! $Id$
65   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE lim_istate
70      !!-------------------------------------------------------------------
71      !!                    ***  ROUTINE lim_istate  ***
72      !!
73      !! ** Purpose :   defined the sea-ice initial state
74      !!
75      !! ** Method  :   
76      !!                This routine will put some ice where ocean
77      !!                is at the freezing point, then fill in ice
78      !!                state variables using prescribed initial
79      !!                values in the namelist           
80      !!
81      !! ** Steps   :   
82      !!                1) Read namelist
83      !!                2) Basal temperature; ice and hemisphere masks
84      !!                3) Fill in the ice thickness distribution using gaussian
85      !!                4) Fill in space-dependent arrays for state variables
86      !!                5) Diagnostic arrays
87      !!                6) Lateral boundary conditions
88      !!
89      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
90      !!              where there is no ice (clem: I do not know why, is it mandatory?)
91      !!
92      !! History :
93      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
94      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
95      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
96      !!--------------------------------------------------------------------
97
98      !! * Local variables
99      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
100      REAL(wp)   :: ztmelts, zdh
101      INTEGER    :: i_hemis, i_fill, jl0 
102      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
103      REAL(wp), POINTER, DIMENSION(:,:)   :: zswitch    ! ice indicator
104      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
105      REAL(wp), POINTER, DIMENSION(:,:)   :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
106      REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini               !data by cattegories to fill
107      !--------------------------------------------------------------------
108
109      CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini )
110      CALL wrk_alloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
111      CALL wrk_alloc( jpi, jpj,      zswitch )
112
113      IF(lwp) WRITE(numout,*)
114      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
115      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
116
117      !--------------------------------------------------------------------
118      ! 1) Read namelist
119      !--------------------------------------------------------------------
120
121      CALL lim_istate_init     !  reading the initials parameters of the ice
122
123      ! surface temperature
124      DO jl = 1, jpl ! loop over categories
125         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
126         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
127      END DO
128
129      ! basal temperature (considered at freezing point)
130      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
131      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
132
133
134      IF( ln_iceini ) THEN
135
136         !--------------------------------------------------------------------
137         ! 2) Basal temperature, ice mask and hemispheric index
138         !--------------------------------------------------------------------
139
140         DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
141            DO ji = 1, jpi
142               IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN
143                  zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
144               ELSE                                                                                   
145                  zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
146               ENDIF
147            END DO
148         END DO
149
150         !--------------------------------------------------------------------
151         ! 3) Initialization of sea ice state variables
152         !--------------------------------------------------------------------
153         IF( ln_iceini_file )THEN
154
155            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
156            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
157            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
158            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
159            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
160            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
161
162         ELSE ! ln_iceini_file = F
163
164            !-----------------------------
165            ! 3.1) Hemisphere-dependent arrays
166            !-----------------------------
167            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
168            DO jj = 1, jpj
169               DO ji = 1, jpi
170                  IF( fcor(ji,jj) >= 0._wp ) THEN
171                     zht_i_ini(ji,jj) = rn_hti_ini_n
172                     zht_s_ini(ji,jj) = rn_hts_ini_n
173                     zat_i_ini(ji,jj) = rn_ati_ini_n
174                     zts_u_ini(ji,jj) = rn_tmi_ini_n
175                     zsm_i_ini(ji,jj) = rn_smi_ini_n
176                     ztm_i_ini(ji,jj) = rn_tmi_ini_n
177                  ELSE
178                     zht_i_ini(ji,jj) = rn_hti_ini_s
179                     zht_s_ini(ji,jj) = rn_hts_ini_s
180                     zat_i_ini(ji,jj) = rn_ati_ini_s
181                     zts_u_ini(ji,jj) = rn_tmi_ini_s
182                     zsm_i_ini(ji,jj) = rn_smi_ini_s
183                     ztm_i_ini(ji,jj) = rn_tmi_ini_s
184                  ENDIF
185               END DO
186            END DO
187
188         ENDIF ! ln_iceini_file
189
190         zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume
191
192         !---------------------------------------------------------------------
193         ! 3.2) Distribute ice concentration and thickness into the categories
194         !---------------------------------------------------------------------
195         ! a gaussian distribution for ice concentration is used
196         ! then we check whether the distribution fullfills
197         ! volume and area conservation, positivity and ice categories bounds
198         zh_i_ini(:,:,:) = 0._wp 
199         za_i_ini(:,:,:) = 0._wp
200         zv_i_ini(:,:,:) = 0._wp
201
202         DO jj = 1, jpj
203            DO ji = 1, jpi
204
205               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
206
207                  ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
208!                  ztests  = 0
209
210                  DO i_fill = jpl, 1, -1
211
212!                     IF( ztests .NE. 4 ) THEN
213                     IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
214                        !----------------------------
215                        ! fill the i_fill categories
216                        !----------------------------
217                        ! *** 1 category to fill
218                        IF ( i_fill .EQ. 1 ) THEN
219                           zh_i_ini(ji,jj,    1)   = zht_i_ini(ji,jj)
220                           za_i_ini(ji,jj,    1)   = zat_i_ini(ji,jj)
221                           zh_i_ini(ji,jj,2:jpl)   = 0._wp
222                           za_i_ini(ji,jj,2:jpl)   = 0._wp
223                        ELSE
224
225                           ! *** >1 categores to fill
226                           !--- Ice thicknesses in the i_fill - 1 first categories
227                           DO jl = 1, i_fill - 1
228                              zh_i_ini(ji,jj,jl) = hi_mean(jl)
229                           END DO
230               
231                           !--- jl0: most likely index where cc will be maximum
232                           DO jl = 1, jpl
233                              IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. &
234                                 & ( zht_i_ini(ji,jj) <= hi_max(jl)   ) ) THEN
235                                 jl0 = jl
236                              ENDIF
237                           END DO
238                           jl0 = MIN(jl0, i_fill)
239               
240                           !--- Concentrations
241                           za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
242                           DO jl = 1, i_fill - 1
243                              IF( jl .NE. jl0 )THEN
244                                 zsigma             = 0.5 * zht_i_ini(ji,jj)
245                                 zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma
246                                 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
247                              ENDIF
248                           END DO
249               
250                           zA = 0. ! sum of the areas in the jpl categories
251                           DO jl = 1, i_fill - 1
252                              zA = zA + za_i_ini(ji,jj,jl)
253                           END DO
254                           za_i_ini(ji,jj,i_fill)   = zat_i_ini(ji,jj) - zA ! ice conc in the last category
255                           IF ( i_fill .LT. jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
256         
257                           !--- Ice thickness in the last category
258                           zV = 0. ! sum of the volumes of the N-1 categories
259                           DO jl = 1, i_fill - 1
260                              zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl)
261                           END DO
262                           zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill) 
263                           IF ( i_fill .LT. jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
264
265                           !--- volumes
266                           zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:)
267                           IF ( i_fill .LT. jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp
268
269                        ENDIF ! i_fill
270
271                        !---------------------
272                        ! Compatibility tests
273                        !---------------------
274                        ! Test 1: area conservation
275                        zA_cons = SUM( za_i_ini(ji,jj,:) )
276                        zconv   = ABS( zat_i_ini(ji,jj) - zA_cons )
277                        IF ( zconv < 1.0e-6 ) THEN
278                           ztest_1 = 1
279                        ELSE 
280                          ! this write is useful
281                          IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,   &
282                           &                                                    ' zat_i_ini = ',zat_i_ini(ji,jj)
283                          ztest_1 = 0
284                        ENDIF
285
286                        ! Test 2: volume conservation
287                        zV_cons = SUM( zv_i_ini(ji,jj,:) )
288                        zconv   = ABS( zvt_i_ini(ji,jj) - zV_cons )
289
290                        IF( zconv < 1.0e-6 ) THEN
291                           ztest_2 = 1
292                        ELSE
293                           ! this write is useful
294                           IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons,   &
295                              &                                                    ' zvt_i_ini = ', zvt_i_ini(ji,jj)
296                           ztest_2 = 0
297                        ENDIF
298
299                        ! Test 3: thickness of the last category is in-bounds ?
300                        IF ( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN
301                           ztest_3 = 1
302                        ELSE
303                           ! this write is useful
304                           IF(lwp) WRITE(numout,*) ' * TEST3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** ',   &
305                              &    ' zh_i_ini(ji,jj,i_fill) = ', zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
306                           IF(lwp) WRITE(numout,*) ' ji,jj,i_fill ',ji,jj,i_fill
307                           IF(lwp) WRITE(numout,*) 'zht_i_ini ',zht_i_ini(ji,jj)
308                           ztest_3 = 0
309                        ENDIF
310
311                        ! Test 4: positivity of ice concentrations
312                        ztest_4 = 1
313                        DO jl = 1, jpl
314                           IF ( za_i_ini(ji,jj,jl) .LT. 0._wp ) THEN 
315                              ! this write is useful
316                              IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(ji,jj,jl)
317                              ztest_4 = 0
318                           ENDIF
319                        END DO
320
321                     ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
322 
323                     ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
324
325                  END DO ! i_fill
326
327                  IF(lwp) THEN
328                     WRITE(numout,*) ' ztests : ', ztests
329                     IF( ztests .NE. 4 )THEN
330                        WRITE(numout,*)
331                        WRITE(numout,*) ' !!!! ALERT                  !!! '
332                        WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
333                        WRITE(numout,*)
334                        WRITE(numout,*) ' *** ztests is not equal to 4 '
335                        WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
336                        WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
337                        WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
338                     ENDIF ! ztests .NE. 4
339                  ENDIF
340     
341               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zhm_i_ini(ji,jj) > 0._wp
342
343            ENDDO   
344         ENDDO   
345
346         !---------------------------------------------------------------------
347         ! 3.3) Space-dependent arrays for ice state variables
348         !---------------------------------------------------------------------
349
350         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
351         DO jl = 1, jpl ! loop over categories
352            DO jj = 1, jpj
353               DO ji = 1, jpi
354                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
355                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
356                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
357                  o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                                    ! age (1 day)
358                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
359
360                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
361                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
362                  ELSE
363                    ht_s(ji,jj,jl)= 0._wp
364                  ENDIF
365
366                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist
367                  ! In case snow load is in excess that would lead to transformation from snow to ice
368                  ! Then, transfer the snow excess into the ice (different from limthd_dh)
369                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
370                  ! recompute ht_i, ht_s avoiding out of bounds values
371                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
372                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
373
374                  ! ice volume, salt content, age content
375                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
376                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
377                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
378                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
379               END DO
380            END DO
381         END DO
382
383         ! Snow temperature and heat content
384         DO jk = 1, nlay_s
385            DO jl = 1, jpl ! loop over categories
386               DO jj = 1, jpj
387                  DO ji = 1, jpi
388                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
389                     ! Snow energy of melting
390                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
391
392                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
393                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
394                  END DO
395               END DO
396            END DO
397         END DO
398
399         ! Ice salinity, temperature and heat content
400         DO jk = 1, nlay_i
401            DO jl = 1, jpl ! loop over categories
402               DO jj = 1, jpj
403                  DO ji = 1, jpi
404                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
405                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
406                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
407
408                     ! heat content per unit volume
409                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
410                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
411                        -   rcp     * ( ztmelts - rt0 ) )
412
413                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
414                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
415                  END DO
416               END DO
417            END DO
418         END DO
419
420         tn_ice (:,:,:) = t_su (:,:,:)
421
422      ELSE ! if ln_iceini=false
423         a_i  (:,:,:) = 0._wp
424         v_i  (:,:,:) = 0._wp
425         v_s  (:,:,:) = 0._wp
426         smv_i(:,:,:) = 0._wp
427         oa_i (:,:,:) = 0._wp
428         ht_i (:,:,:) = 0._wp
429         ht_s (:,:,:) = 0._wp
430         sm_i (:,:,:) = 0._wp
431         o_i  (:,:,:) = 0._wp
432
433         e_i(:,:,:,:) = 0._wp
434         e_s(:,:,:,:) = 0._wp
435
436         DO jl = 1, jpl
437            DO jk = 1, nlay_i
438               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
439            END DO
440            DO jk = 1, nlay_s
441               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
442            END DO
443         END DO
444
445      ENDIF ! ln_iceini
446     
447      at_i (:,:) = 0.0_wp
448      DO jl = 1, jpl
449         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
450      END DO
451      !
452      !--------------------------------------------------------------------
453      ! 4) Global ice variables for output diagnostics                    |
454      !--------------------------------------------------------------------
455      u_ice (:,:)     = 0._wp
456      v_ice (:,:)     = 0._wp
457      stress1_i(:,:)  = 0._wp
458      stress2_i(:,:)  = 0._wp
459      stress12_i(:,:) = 0._wp
460
461      !--------------------------------------------------------------------
462      ! 5) Moments for advection
463      !--------------------------------------------------------------------
464
465      sxopw (:,:) = 0._wp 
466      syopw (:,:) = 0._wp 
467      sxxopw(:,:) = 0._wp 
468      syyopw(:,:) = 0._wp 
469      sxyopw(:,:) = 0._wp
470
471      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
472      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
473      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
474      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
475      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
476
477      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
478      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
479      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
480      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
481      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
482
483      sxsal  (:,:,:)  = 0._wp
484      sysal  (:,:,:)  = 0._wp
485      sxxsal (:,:,:)  = 0._wp
486      syysal (:,:,:)  = 0._wp
487      sxysal (:,:,:)  = 0._wp
488
489      sxage  (:,:,:)  = 0._wp
490      syage  (:,:,:)  = 0._wp
491      sxxage (:,:,:)  = 0._wp
492      syyage (:,:,:)  = 0._wp
493      sxyage (:,:,:)  = 0._wp
494
495
496      CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini,  za_i_ini,  zv_i_ini )
497      CALL wrk_dealloc( jpi, jpj,      zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
498      CALL wrk_dealloc( jpi, jpj,      zswitch )
499
500   END SUBROUTINE lim_istate
501
502   SUBROUTINE lim_istate_init
503      !!-------------------------------------------------------------------
504      !!                   ***  ROUTINE lim_istate_init  ***
505      !!       
506      !! ** Purpose : Definition of initial state of the ice
507      !!
508      !! ** Method : Read the namiceini namelist and check the parameter
509      !!       values called at the first timestep (nit000)
510      !!
511      !! ** input :
512      !!        Namelist namiceini
513      !!
514      !! history :
515      !!  8.5  ! 03-08 (C. Ethe) original code
516      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
517      !!-----------------------------------------------------------------------------
518      !
519      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
520      INTEGER :: ji,jj
521      INTEGER :: ifpr, ierror
522      !
523      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ice files
524      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
525      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
526      !
527      NAMELIST/namiceini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
528         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
529         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
530         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
531      !!-----------------------------------------------------------------------------
532      !
533      REWIND( numnam_ice_ref )         ! Namelist namiceini in reference namelist : Ice initial state
534      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
535901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
536
537      REWIND( numnam_ice_cfg )         ! Namelist namiceini in configuration namelist : Ice initial state
538      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
539902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
540      IF(lwm) WRITE ( numoni, namiceini )
541
542      slf_i(jp_hti) = sn_hti   ;   slf_i(jp_hts) = sn_hts
543      slf_i(jp_ati) = sn_ati   ;   slf_i(jp_tsu) = sn_tsu
544      slf_i(jp_tmi) = sn_tmi   ;   slf_i(jp_smi) = sn_smi
545
546      IF(lwp) THEN                     ! control print
547         WRITE(numout,*)
548         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
549         WRITE(numout,*) '~~~~~~~~~~~~~~~'
550         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
551         WRITE(numout,*) '   ice initialization from a netcdf file      ln_iceini_file  = ', ln_iceini_file
552         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
553         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
554         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
555         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
556         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
557         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
558         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
559         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
560         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
561         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
562         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
563      ENDIF
564
565      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
566         !
567         ! set si structure
568         ALLOCATE( si(jpfldi), STAT=ierror )
569         IF( ierror > 0 ) THEN
570            CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' )   ;   RETURN
571         ENDIF
572         !
573         DO ifpr = 1, jpfldi
574            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
575            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
576         END DO
577         !
578         ! fill si with slf_i and control print
579         CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' )
580         !
581         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
582         !
583      ENDIF
584      !
585   END SUBROUTINE lim_istate_init
586
587#else
588   !!----------------------------------------------------------------------
589   !!   Default option :         Empty module          NO LIM sea-ice model
590   !!----------------------------------------------------------------------
591CONTAINS
592   SUBROUTINE lim_istate          ! Empty routine
593   END SUBROUTINE lim_istate
594#endif
595
596   !!======================================================================
597END MODULE limistate
Note: See TracBrowser for help on using the repository browser.