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/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/iceistate.F90 @ 12489

Last change on this file since 12489 was 12489, checked in by davestorkey, 4 years ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp -> rn_atfp (namelist parameter used everywhere) and rau0 -> rho0.

This version gives bit-comparable results to the previous version of the trunk.

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