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.
iceistate.F90 in NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE – NEMO

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/iceistate.F90 @ 11402

Last change on this file since 11402 was 11402, checked in by clem, 5 years ago

debug previous commit
  • Property svn:keywords set to Id
File size: 25.8 KB
Line 
1MODULE iceistate
2   !!======================================================================
3   !!                     ***  MODULE  iceistate  ***
4   !!   sea-ice : Initialization of ice variables
5   !!======================================================================
6   !! History :  2.0  !  2004-01  (C. Ethe, G. Madec) Original code
7   !!            3.0  !  2007     (M. Vancoppenolle)  Rewrite for ice cats
8   !!            4.0  !  2018     (many people)       SI3 [aka Sea Ice cube]
9   !!----------------------------------------------------------------------
10#if defined key_si3
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                       SI3 sea-ice model
13   !!----------------------------------------------------------------------
14   !!   ice_istate       :  initialization of diagnostics ice variables
15   !!   ice_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 , ONLY : sst_m, sss_m, ln_ice_embd 
21   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
22   USE eosbn2         ! equation of state
23   USE domvvl         ! Variable volume
24   USE ice            ! sea-ice: variables
25   USE ice1D          ! sea-ice: thermodynamics variables
26   USE icetab         ! sea-ice: 1D <==> 2D transformation
27   USE icevar         ! sea-ice: operations
28   !
29   USE in_out_manager ! I/O manager
30   USE iom            ! I/O manager library
31   USE lib_mpp        ! MPP library
32   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
33   USE fldread        ! read input fields
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   ice_istate        ! called by icestp.F90
39   PUBLIC   ice_istate_init   ! called by icestp.F90
40   !
41   !                             !! ** namelist (namini) **
42   LOGICAL, PUBLIC  ::   ln_iceini        !: Ice initialization or not
43   LOGICAL, PUBLIC  ::   ln_iceini_file   !: Ice initialization from 2D netcdf file
44   REAL(wp) ::   rn_thres_sst
45   REAL(wp) ::   rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
46   REAL(wp) ::   rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
47   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n
48   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s
49   !
50   !                              ! if ln_iceini_file = T
51   INTEGER , PARAMETER ::   jpfldi = 9           ! maximum number of files to read
52   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m)
53   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m)
54   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction     (-)
55   INTEGER , PARAMETER ::   jp_smi = 4           ! index of ice salinity     (g/kg)
56   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temperature  (K)
57   INTEGER , PARAMETER ::   jp_tsu = 6           ! index of ice surface temp (K)
58   INTEGER , PARAMETER ::   jp_tms = 7           ! index of snw temperature  (K)
59   INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-)
60   INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
62   !   
63   !!----------------------------------------------------------------------
64   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
65   !! $Id$
66   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE ice_istate( kt )
71      !!-------------------------------------------------------------------
72      !!                    ***  ROUTINE ice_istate  ***
73      !!
74      !! ** Purpose :   defined the sea-ice initial state
75      !!
76      !! ** Method  :   This routine will put some ice where ocean
77      !!                is at the freezing point, then fill in ice
78      !!                state variables using prescribed initial
79      !!                values in the namelist           
80      !!
81      !! ** Steps   :   1) Set initial surface and basal temperatures
82      !!                2) Recompute or read sea ice state variables
83      !!                3) Fill in the ice thickness distribution using gaussian
84      !!                4) Fill in space-dependent arrays for state variables
85      !!                5) snow-ice mass computation
86      !!                6) store before fields
87      !!
88      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
89      !!              where there is no ice
90      !!--------------------------------------------------------------------
91      INTEGER, INTENT(in) ::   kt   ! time step
92      !!
93      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices
94      REAL(wp) ::   ztmelts
95      INTEGER , DIMENSION(4)           ::   itest
96      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
97      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
98      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file
99      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
100      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file
101      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays
102      !!
103      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d
104      !--------------------------------------------------------------------
105
106      IF(lwp) WRITE(numout,*)
107      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
108      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
109
110      !---------------------------
111      ! 1) 1st init. of the fields
112      !---------------------------
113      !
114      ! basal temperature (considered at freezing point)   [Kelvin]
115      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
116      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
117      !
118      ! surface temperature and conductivity
119      DO jl = 1, jpl
120         t_su   (:,:,jl) = rt0 * tmask(:,:,1)  ! temp at the surface
121         cnd_ice(:,:,jl) = 0._wp               ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
122      END DO
123      !
124      ! ice and snw temperatures
125      DO jl = 1, jpl
126         DO jk = 1, nlay_i
127            t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
128         END DO
129         DO jk = 1, nlay_s
130            t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
131         END DO
132      END DO
133      !
134      ! specific temperatures for coupled runs
135      tn_ice (:,:,:) = t_i (:,:,1,:)
136      t1_ice (:,:,:) = t_i (:,:,1,:)
137
138      ! heat contents
139      e_i (:,:,:,:) = 0._wp
140      e_s (:,:,:,:) = 0._wp
141     
142      ! general fields
143      a_i (:,:,:) = 0._wp
144      v_i (:,:,:) = 0._wp
145      v_s (:,:,:) = 0._wp
146      sv_i(:,:,:) = 0._wp
147      oa_i(:,:,:) = 0._wp
148      !
149      h_i (:,:,:) = 0._wp
150      h_s (:,:,:) = 0._wp
151      s_i (:,:,:) = 0._wp
152      o_i (:,:,:) = 0._wp
153      !
154      ! melt ponds
155      a_ip     (:,:,:) = 0._wp
156      v_ip     (:,:,:) = 0._wp
157      a_ip_frac(:,:,:) = 0._wp
158      h_ip     (:,:,:) = 0._wp
159      !
160      ! ice velocities
161      u_ice (:,:) = 0._wp
162      v_ice (:,:) = 0._wp
163      !
164      !------------------------------------------------------------------------
165      ! 2) overwrite some of the fields with namelist parameters or netcdf file
166      !------------------------------------------------------------------------
167      IF( ln_iceini ) THEN
168         !                             !---------------!
169         IF( ln_iceini_file )THEN      ! Read a file   !
170            !                          !---------------!
171            WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp
172            ELSEWHERE                     ;   zswitch(:,:) = 0._wp
173            END WHERE
174            !
175            CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
176            !
177            ! -- mandatory fields -- !
178            zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
179            zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
180            zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
181
182            ! -- optional fields -- !
183            !    if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
184            !
185            ! ice salinity
186            IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
187               &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
188            zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
189            !
190            ! ice temperature
191            IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) &
192               &     si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
193            ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
194            !
195            ! surface temperature => set to ice temperature if it exists
196            IF    ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
197                     si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
198            ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
199                     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
200            ENDIF
201            zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
202            !
203            ! snow temperature => set to ice temperature if it exists
204            IF    ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
205                     si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
206            ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
207                     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
208            ENDIF
209            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1)
210            !
211            ! ponds
212            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
213               &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
214            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
215            IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
216               &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
217            zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
218            !
219            ! change the switch for the following
220            WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
221            ELSEWHERE                         ;   zswitch(:,:) = 0._wp
222            END WHERE
223            !                          !---------------!
224         ELSE                          ! Read namelist !
225            !                          !---------------!
226            ! no ice if (sst - Tfreez) >= thresold
227            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
228            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
229            END WHERE
230            !
231            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
232            WHERE( ff_t(:,:) >= 0._wp )
233               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
234               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
235               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
236               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
237               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
238               zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
239               ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
240               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:)
241               zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
242            ELSEWHERE
243               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
244               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
245               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
246               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
247               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
248               zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
249               ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
250               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:)
251               zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
252            END WHERE
253            !
254         ENDIF
255         !-------------!
256         ! fill fields !
257         !-------------!
258         ! select ice covered grid points
259         npti = 0 ; nptidx(:) = 0
260         DO jj = 1, jpj
261            DO ji = 1, jpi
262               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
263                  npti         = npti  + 1
264                  nptidx(npti) = (jj - 1) * jpi + ji
265               ENDIF
266            END DO
267         END DO
268
269         ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
270         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini )
271         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini )
272         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini )
273         CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
274         CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
275         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini )
276         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini )
277
278         ! allocate temporary arrays
279         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
280            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl) )
281         
282         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
283         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                   zhi_2d, zhs_2d, zai_2d , &
284            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti),   zti_2d, zts_2d, ztsu_2d, zsi_2d )
285
286         ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
287         DO jl = 1, jpl
288            zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
289            zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
290         END DO
291         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d  , h_i  )
292         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d  , h_s  )
293         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d  , a_i  )
294         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d  , zti_3d )
295         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d  , zts_3d )
296         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
297         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d  , s_i  )
298
299         ! deallocate temporary arrays
300         DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
301            &        zti_2d, zts_2d, ztsu_2d, zsi_2d )
302
303         ! Melt ponds: distribute uniformely over the categories
304         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN
305            DO jl = 1, jpl
306               a_ip_frac(:,:,jl) = zapnd_ini(:,:)
307               h_ip     (:,:,jl) = zhpnd_ini(:,:)
308               a_ip     (:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
309               v_ip     (:,:,jl) = h_ip     (:,:,jl) * a_ip(:,:,jl)
310            END DO
311         ENDIF
312         
313         ! calculate extensive and intensive variables
314         CALL ice_var_salprof ! for sz_i
315         DO jl = 1, jpl
316            DO jj = 1, jpj
317               DO ji = 1, jpi
318                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
319                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
320                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
321               END DO
322            END DO
323         END DO
324         !
325         DO jl = 1, jpl
326            DO jk = 1, nlay_s
327               DO jj = 1, jpj
328                  DO ji = 1, jpi
329                     t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
330                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
331                        &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
332                  END DO
333               END DO
334            END DO
335         END DO
336         !
337         DO jl = 1, jpl
338            DO jk = 1, nlay_i
339               DO jj = 1, jpj
340                  DO ji = 1, jpi
341                     t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
342                     ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
343                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
344                        &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
345                        &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
346                        &                       - rcp   * ( ztmelts - rt0 ) )
347                  END DO
348               END DO
349            END DO
350         END DO
351
352         ! specific temperatures for coupled runs
353         tn_ice(:,:,:) = t_su(:,:,:)
354         t1_ice(:,:,:) = t_i (:,:,1,:)
355         !
356      ENDIF ! ln_iceini
357      !
358      at_i(:,:) = SUM( a_i, dim=3 )
359      !
360      !----------------------------------------------
361      ! 3) Snow-ice mass (case ice is fully embedded)
362      !----------------------------------------------
363      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass
364      snwice_mass_b(:,:) = snwice_mass(:,:)
365      !
366      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
367         !
368         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
369         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
370         !
371         IF( .NOT.ln_linssh ) THEN
372            !
373            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
374            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
375            !
376            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors               
377               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
378               e3t_b(:,:,jk) = e3t_n(:,:,jk)
379               e3t_a(:,:,jk) = e3t_n(:,:,jk)
380            END DO
381            !
382            ! Reconstruction of all vertical scale factors at now and before time-steps
383            ! =========================================================================
384            ! Horizontal scale factor interpolations
385            ! --------------------------------------
386            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
387            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
388            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
389            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
390            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
391            ! Vertical scale factor interpolations
392            ! ------------------------------------
393            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
394            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
395            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
396            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
397            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
398            ! t- and w- points depth
399            ! ----------------------
400            !!gm not sure of that....
401            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
402            gdepw_n(:,:,1) = 0.0_wp
403            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
404            DO jk = 2, jpk
405               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
406               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
407               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
408            END DO
409         ENDIF
410      ENDIF
411     
412      !------------------------------------
413      ! 4) store fields at before time-step
414      !------------------------------------
415      ! it is only necessary for the 1st interpolation by Agrif
416      a_i_b  (:,:,:)   = a_i  (:,:,:)
417      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
418      v_i_b  (:,:,:)   = v_i  (:,:,:)
419      v_s_b  (:,:,:)   = v_s  (:,:,:)
420      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
421      sv_i_b (:,:,:)   = sv_i (:,:,:)
422      oa_i_b (:,:,:)   = oa_i (:,:,:)
423      u_ice_b(:,:)     = u_ice(:,:)
424      v_ice_b(:,:)     = v_ice(:,:)
425      ! total concentration is needed for Lupkes parameterizations
426      at_i_b (:,:)     = at_i (:,:) 
427
428!!clem: output of initial state should be written here but it is impossible because
429!!      the ocean and ice are in the same file
430!!      CALL dia_wri_state( 'output.init' )
431      !
432   END SUBROUTINE ice_istate
433
434
435   SUBROUTINE ice_istate_init
436      !!-------------------------------------------------------------------
437      !!                   ***  ROUTINE ice_istate_init  ***
438      !!       
439      !! ** Purpose :   Definition of initial state of the ice
440      !!
441      !! ** Method  :   Read the namini namelist and check the parameter
442      !!              values called at the first timestep (nit000)
443      !!
444      !! ** input   :  Namelist namini
445      !!
446      !!-----------------------------------------------------------------------------
447      INTEGER ::   ios   ! Local integer output status for namelist read
448      INTEGER ::   ifpr, ierror
449      !
450      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
451      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
452      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
453      !
454      NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
455         &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
456         &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
457         &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
458         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
459         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
460      !!-----------------------------------------------------------------------------
461      !
462      REWIND( numnam_ice_ref )              ! Namelist namini in reference namelist : Ice initial state
463      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
464901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namini in reference namelist' )
465      REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state
466      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
467902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namini in configuration namelist' )
468      IF(lwm) WRITE ( numoni, namini )
469      !
470      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
471      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi
472      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms
473      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd
474      !
475      IF(lwp) THEN                          ! control print
476         WRITE(numout,*)
477         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
478         WRITE(numout,*) '~~~~~~~~~~~~~~~'
479         WRITE(numout,*) '   Namelist namini:'
480         WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini
481         WRITE(numout,*) '      ice initialization from a netcdf file            ln_iceini_file = ', ln_iceini_file
482         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst
483         IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
484            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
485            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s
486            WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s
487            WRITE(numout,*) '      initial ice salinity  in the north-south         rn_smi_ini     = ', rn_smi_ini_n,rn_smi_ini_s
488            WRITE(numout,*) '      initial surf temperat in the north-south         rn_tsu_ini     = ', rn_tsu_ini_n,rn_tsu_ini_s
489            WRITE(numout,*) '      initial ice temperat  in the north-south         rn_tmi_ini     = ', rn_tmi_ini_n,rn_tmi_ini_s
490            WRITE(numout,*) '      initial snw temperat  in the north-south         rn_tms_ini     = ', rn_tms_ini_n,rn_tms_ini_s
491            WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s
492            WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s
493         ENDIF
494      ENDIF
495      !
496      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
497         !
498         ! set si structure
499         ALLOCATE( si(jpfldi), STAT=ierror )
500         IF( ierror > 0 ) THEN
501            CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' )   ;   RETURN
502         ENDIF
503         !
504         DO ifpr = 1, jpfldi
505            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
506            IF( slf_i(ifpr)%ln_tint )  ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
507         END DO
508         !
509         ! fill si with slf_i and control print
510         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
511         !
512      ENDIF
513      !
514   END SUBROUTINE ice_istate_init
515
516#else
517   !!----------------------------------------------------------------------
518   !!   Default option :         Empty module         NO SI3 sea-ice model
519   !!----------------------------------------------------------------------
520#endif
521
522   !!======================================================================
523END MODULE iceistate
Note: See TracBrowser for help on using the repository browser.