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

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 @ 8061

Last change on this file since 8061 was 8061, checked in by vancop, 7 years ago

Quick commit on empirical melt ponds

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