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

source: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 6723

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

#1751 - branch SIMPLIF_6_aerobulk: add aerobulk package including NCAR, COARE and ECMWF bulk

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