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/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_retain_melt/src/ICE/iceistate.F90 @ 13080

Last change on this file since 13080 was 13080, checked in by dancopsey, 4 years ago

Merge in existing fix_cpl branch.

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