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/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceistate.F90 @ 13350

Last change on this file since 13350 was 13284, checked in by smasson, 4 years ago
4.0-HEAD: merge 4.0-HEAD_r12713_clem_dan_fixcpl into 4.0-HEAD
Property svn:keywords set to `Id`
File size: 28.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: