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

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 7917

Last change on this file since 7917 was 7917, checked in by timgraham, 7 years ago

A few corrections for SETTE compilation

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