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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/iceistate.F90 @ 11822

Last change on this file since 11822 was 11822, checked in by acc, 4 years ago

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