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_r6522_SIMPLIF_3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2016/dev_r6522_SIMPLIF_3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 6863

Last change on this file since 6863 was 6470, checked in by clem, 8 years ago

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