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/2020/SI3-03_VP_rheology/src/ICE – NEMO

source: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/iceistate.F90 @ 13681

Last change on this file since 13681 was 13681, checked in by vancop, 3 years ago

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