Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

sbcice_cice.F90

Last change on this file was 9828, checked in by dancopsey, 6 years ago
Stop multiplying zevap_ice by -1.
File size: 46.8 KB

Rev	Line
[2874]	1	MODULE sbcice_cice
	2	!!======================================================================
	3	!! * MODULE sbcice_cice *
	4	!! To couple with sea ice model CICE (LANL)
	5	!!=====================================================================
	6	#if defined key_cice
	7	!!----------------------------------------------------------------------
	8	!! 'key_cice' : CICE sea-ice model
	9	!!----------------------------------------------------------------------
	10	!! sbc_ice_cice : sea-ice model time-stepping and update ocean sbc over ice-covered area
	11	!!----------------------------------------------------------------------
	12	USE oce ! ocean dynamics and tracers
	13	USE dom_oce ! ocean space and time domain
[3275]	14	USE domvvl
[9500]	15	USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater
	16	USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0
[2874]	17	USE in_out_manager ! I/O manager
[4990]	18	USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit
[2874]	19	USE lib_mpp ! distributed memory computing library
	20	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	21	USE daymod ! calendar
	22	USE fldread ! read input fields
	23	USE sbc_oce ! Surface boundary condition: ocean fields
	24	USE sbc_ice ! Surface boundary condition: ice fields
[7646]	25	USE sbcblk ! Surface boundary condition: bulk
[2874]	26	USE sbccpl
	27
	28	USE ice_kinds_mod
	29	USE ice_blocks
	30	USE ice_domain
	31	USE ice_domain_size
	32	USE ice_boundary
	33	USE ice_constants
	34	USE ice_gather_scatter
	35	USE ice_calendar, only: dt
[9500]	36	# if defined key_cice4
[9498]	37	USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
[2874]	38	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[5133]	39	strocnxT,strocnyT, &
[3189]	40	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, &
	41	fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, &
[2874]	42	flatn_f,fsurfn_f,fcondtopn_f, &
	43	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
[9500]	44	swvdr,swvdf,swidr,swidf,Tf
[4990]	45	USE ice_therm_vertical, only: calc_Tsfc
	46	#else
[9500]	47	USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,&
	48	vsnon,vice,vicen,nt_Tsfc
[4990]	49	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[5133]	50	strocnxT,strocnyT, &
[9500]	51	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai, &
	52	fresh_ai,fhocn_ai,fswthru_ai,frzmlt, &
[4990]	53	flatn_f,fsurfn_f,fcondtopn_f, &
[9500]	54	#ifdef key_asminc
	55	daice_da,fresh_da,fsalt_da, &
	56	#endif
[4990]	57	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
[9500]	58	swvdr,swvdf,swidr,swidf,Tf, &
	59	!! When using NEMO with CICE, this change requires use of
	60	!! one of the following two CICE branches:
	61	!! - at CICE5.0, hadax/r1015_GSI8_with_GSI7
	62	!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
	63	keffn_top,Tn_top
	64
	65	USE ice_therm_shared, only: calc_Tsfc, heat_capacity
	66	USE ice_shortwave, only: apeffn
[4990]	67	#endif
[2874]	68	USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
[3176]	69	USE ice_atmo, only: calc_strair
[2874]	70
	71	USE CICE_InitMod
	72	USE CICE_RunMod
	73	USE CICE_FinalMod
	74
	75	IMPLICIT NONE
	76	PRIVATE
	77
	78	PUBLIC cice_sbc_init ! routine called by sbc_init
	79	PUBLIC cice_sbc_final ! routine called by sbc_final
	80	PUBLIC sbc_ice_cice ! routine called by sbc
	81
[4627]	82	INTEGER :: ji_off
	83	INTEGER :: jj_off
[3625]	84
[2874]	85	INTEGER , PARAMETER :: jpfld = 13 ! maximum number of files to read
	86	INTEGER , PARAMETER :: jp_snow = 1 ! index of snow file
	87	INTEGER , PARAMETER :: jp_rain = 2 ! index of rain file
	88	INTEGER , PARAMETER :: jp_sblm = 3 ! index of sublimation file
	89	INTEGER , PARAMETER :: jp_top1 = 4 ! index of category 1 topmelt file
	90	INTEGER , PARAMETER :: jp_top2 = 5 ! index of category 2 topmelt file
	91	INTEGER , PARAMETER :: jp_top3 = 6 ! index of category 3 topmelt file
	92	INTEGER , PARAMETER :: jp_top4 = 7 ! index of category 4 topmelt file
	93	INTEGER , PARAMETER :: jp_top5 = 8 ! index of category 5 topmelt file
	94	INTEGER , PARAMETER :: jp_bot1 = 9 ! index of category 1 botmelt file
	95	INTEGER , PARAMETER :: jp_bot2 = 10 ! index of category 2 botmelt file
	96	INTEGER , PARAMETER :: jp_bot3 = 11 ! index of category 3 botmelt file
	97	INTEGER , PARAMETER :: jp_bot4 = 12 ! index of category 4 botmelt file
	98	INTEGER , PARAMETER :: jp_bot5 = 13 ! index of category 5 botmelt file
	99	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
	100
	101	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE :: png ! local array used in sbc_cice_ice
	102
[5836]	103	!!----------------------------------------------------------------------
	104	!! NEMO/OPA 3.7 , NEMO-consortium (2015)
[5215]	105	!! $Id$
[5836]	106	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	107	!!----------------------------------------------------------------------
[2874]	108	CONTAINS
	109
	110	INTEGER FUNCTION sbc_ice_cice_alloc()
	111	!!----------------------------------------------------------------------
	112	!! * FUNCTION sbc_ice_cice_alloc *
	113	!!----------------------------------------------------------------------
	114	ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
	115	IF( lk_mpp ) CALL mpp_sum ( sbc_ice_cice_alloc )
	116	IF( sbc_ice_cice_alloc > 0 ) CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
	117	END FUNCTION sbc_ice_cice_alloc
	118
[4990]	119	SUBROUTINE sbc_ice_cice( kt, ksbc )
[2874]	120	!!---------------------------------------------------------------------
	121	!! * ROUTINE sbc_ice_cice *
	122	!!
	123	!! ** Purpose : update the ocean surface boundary condition via the
	124	!! CICE Sea Ice Model time stepping
	125	!!
[3040]	126	!! ** Method : - Get any extra forcing fields for CICE
	127	!! - Prepare forcing fields
[2874]	128	!! - CICE model time stepping
	129	!! - call the routine that computes mass and
	130	!! heat fluxes at the ice/ocean interface
	131	!!
	132	!! ** Action : - time evolution of the CICE sea-ice model
	133	!! - update all sbc variables below sea-ice:
[3625]	134	!! utau, vtau, qns , qsr, emp , sfx
[2874]	135	!!---------------------------------------------------------------------
	136	INTEGER, INTENT(in) :: kt ! ocean time step
[4990]	137	INTEGER, INTENT(in) :: ksbc ! surface forcing type
[2874]	138	!!----------------------------------------------------------------------
[3193]	139	!
[2874]	140	! !----------------------!
	141	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Ice time-step only !
	142	! !----------------------!
	143
	144	! Make sure any fluxes required for CICE are set
[4990]	145	IF ( ksbc == jp_flx ) THEN
[2874]	146	CALL cice_sbc_force(kt)
[5407]	147	ELSE IF ( ksbc == jp_purecpl ) THEN
[9019]	148	CALL sbc_cpl_ice_flx( fr_i )
[2874]	149	ENDIF
	150
[4990]	151	CALL cice_sbc_in ( kt, ksbc )
[2874]	152	CALL CICE_Run
[4990]	153	CALL cice_sbc_out ( kt, ksbc )
[2874]	154
[5407]	155	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(kt+1)
[2874]	156
	157	ENDIF ! End sea-ice time step only
[3193]	158	!
[2874]	159	END SUBROUTINE sbc_ice_cice
	160
[5836]	161
	162	SUBROUTINE cice_sbc_init( ksbc )
[2874]	163	!!---------------------------------------------------------------------
	164	!! * ROUTINE cice_sbc_init *
[3040]	165	!! ** Purpose: Initialise ice related fields for NEMO and coupling
[2874]	166	!!
[5836]	167	!!---------------------------------------------------------------------
[4990]	168	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
[9125]	169	REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
[9500]	170	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d
[4990]	171	INTEGER :: ji, jj, jl, jk ! dummy loop indices
[2874]	172	!!---------------------------------------------------------------------
[3193]	173	!
[2874]	174	IF(lwp) WRITE(numout,*)'cice_sbc_init'
	175
[4627]	176	ji_off = INT ( (jpiglo - nx_global) / 2 )
	177	jj_off = INT ( (jpjglo - ny_global) / 2 )
	178
[9500]	179	! Initialize CICE
[3176]	180	CALL CICE_Initialize
[2874]	181
[9500]	182	! Do some CICE consistency checks
[5407]	183	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
[3193]	184	IF ( calc_strair .OR. calc_Tsfc ) THEN
	185	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
	186	ENDIF
[7646]	187	ELSEIF (ksbc == jp_blk) THEN
[3193]	188	IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
	189	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
	190	ENDIF
	191	ENDIF
[3176]	192
	193
[9500]	194	! allocate sbc_ice and sbc_cice arrays
	195	IF( sbc_ice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_alloc : unable to allocate arrays' )
[2874]	196	IF( sbc_ice_cice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
	197
[9500]	198	! Ensure that no temperature points are below freezing if not a NEMO restart
[2874]	199	IF( .NOT. ln_rstart ) THEN
[9500]	200
	201	DO jk=1,jpk
	202	CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk) )
	203	ENDDO
	204	tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d )
[2874]	205	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
	206
[9500]	207	#if defined key_nemocice_decomp
	208	! Pass initial SST from NEMO to CICE so ice is initialised correctly if
	209	! there is no restart file.
	210	! Values from a CICE restart file would overwrite this
[9819]	211	CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.,'sst1')
[9500]	212	#endif
[2874]	213
[9500]	214	ENDIF
	215
	216	! calculate surface freezing temperature and send to CICE
	217	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1))
[9819]	218	CALL nemo2cice(sstfrz,Tf, 'T', 1.,'tf1' )
[9500]	219
[3193]	220	CALL cice2nemo(aice,fr_i, 'T', 1. )
[5407]	221	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
[3625]	222	DO jl=1,ncat
	223	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	224	ENDDO
	225	ENDIF
[2874]	226
	227	! T point to U point
	228	! T point to V point
[9500]	229	fr_iu(:,:)=0.0
	230	fr_iv(:,:)=0.0
[3193]	231	DO jj=1,jpjm1
	232	DO ji=1,jpim1
	233	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	234	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	235	ENDDO
	236	ENDDO
[2874]	237
[9098]	238	CALL lbc_lnk_multi( fr_iu , 'U', 1., fr_iv , 'V', 1. )
[3625]	239
[9019]	240	! set the snow+ice mass
	241	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	242	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	243	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	244	snwice_mass_b(:,:) = snwice_mass(:,:)
	245
[6140]	246	IF( .NOT.ln_rstart ) THEN
[9019]	247	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
[4990]	248	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	249	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
[6140]	250
	251	!!gm This should be put elsewhere.... (same remark for limsbc)
	252	!!gm especially here it is assumed zstar coordinate, but it can be ztilde....
	253	IF( .NOT.ln_linssh ) THEN
	254	!
	255	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	256	e3t_n(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshn(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
	257	e3t_b(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshb(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
	258	ENDDO
	259	e3t_a(:,:,:) = e3t_b(:,:,:)
	260	! Reconstruction of all vertical scale factors at now and before time-steps
	261	! =============================================================================
	262	! Horizontal scale factor interpolations
	263	! --------------------------------------
	264	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
	265	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
	266	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	267	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
	268	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
	269	! Vertical scale factor interpolations
	270	! ------------------------------------
	271	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
	272	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	273	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
	274	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	275	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
	276	! t- and w- points depth
	277	! ----------------------
	278	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	279	gdepw_n(:,:,1) = 0.0_wp
	280	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
	281	DO jk = 2, jpk
	282	gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk)
	283	gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
	284	gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:)
	285	END DO
	286	ENDIF
[4990]	287	ENDIF
[3625]	288	ENDIF
[9500]	289
	290	#if defined key_asminc
	291	! Initialize fresh water and salt fluxes from data assim
	292	! and data assimilation index to cice
	293	nfresh_da(:,:) = 0.0
	294	nfsalt_da(:,:) = 0.0
	295	ndaice_da(:,:) = 0.0
	296	#endif
[6140]	297	!
[9500]	298	! In coupled mode get extra fields from CICE for passing back to atmosphere
	299
	300	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(nit000)
	301	!
[2874]	302	END SUBROUTINE cice_sbc_init
	303
[3152]	304
[5836]	305	SUBROUTINE cice_sbc_in( kt, ksbc )
[2874]	306	!!---------------------------------------------------------------------
	307	!! * ROUTINE cice_sbc_in *
[3040]	308	!! ** Purpose: Set coupling fields and pass to CICE
[2874]	309	!!---------------------------------------------------------------------
[3152]	310	INTEGER, INTENT(in ) :: kt ! ocean time step
[4990]	311	INTEGER, INTENT(in ) :: ksbc ! surface forcing type
[5836]	312	!
[3625]	313	INTEGER :: ji, jj, jl ! dummy loop indices
[9125]	314	REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zpice
	315	REAL(wp), DIMENSION(jpi,jpj,ncat) :: ztmpn
[3625]	316	REAL(wp) :: zintb, zintn ! dummy argument
[3152]	317	!!---------------------------------------------------------------------
[3193]	318	!
	319	IF( kt == nit000 ) THEN
[2874]	320	IF(lwp) WRITE(numout,*)'cice_sbc_in'
[3193]	321	ENDIF
[2874]	322
[3193]	323	ztmp(:,:)=0.0
[2874]	324
	325	! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
	326	! the first time-step)
	327
	328	! forced and coupled case
	329
[5407]	330	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
[2874]	331
[3193]	332	ztmpn(:,:,:)=0.0
[2874]	333
	334	! x comp of wind stress (CI_1)
	335	! U point to F point
[3193]	336	DO jj=1,jpjm1
	337	DO ji=1,jpi
	338	ztmp(ji,jj) = 0.5 * ( fr_iu(ji,jj) * utau(ji,jj) &
	339	+ fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
	340	ENDDO
	341	ENDDO
[9819]	342
	343	CALL nemo2cice(ztmp,strax,'F', -1.,'strax' )
[2874]	344
	345	! y comp of wind stress (CI_2)
	346	! V point to F point
[3193]	347	DO jj=1,jpj
	348	DO ji=1,jpim1
	349	ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
	350	+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
	351	ENDDO
	352	ENDDO
[9819]	353	CALL nemo2cice(ztmp,stray,'F', -1.,'stray' )
[2874]	354
[9500]	355
	356	! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in
	357	! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby
	358	! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing
	359	! gridbox mean fluxes in the UM by future ice concentration obtained through
	360	! OASIS. This allows for a much more realistic apportionment of energy through
	361	! the ice - and conserves energy.
	362	! Therefore the fluxes are now divided by ice concentration in the coupled
	363	! formulation (jp_purecpl) as well as for jp_flx. This NEMO branch should only
	364	! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at
	365	! which point the GSI8 UM changes were committed.
	366
[2874]	367	! Surface downward latent heat flux (CI_5)
[4990]	368	IF (ksbc == jp_flx) THEN
[3625]	369	DO jl=1,ncat
[9828]	370	ztmpn(:,:,jl)=zevap_ice(:,:,1)a_i(:,:,jl) Lsub
[3193]	371	ENDDO
[9500]	372	ELSE IF (ksbc == jp_purecpl) THEN
	373	DO jl=1,ncat
[9828]	374	ztmpn(:,:,jl)=zevap_ice(:,:,jl)a_i(:,:,jl) Lsub
[3193]	375	ENDDO
[9500]	376	ELSE
[9675]	377	!In coupled mode - zevap_ice calculated in sbc_cpl for each category
[9828]	378	ztmpn(:,:,1:ncat)=zevap_ice(:,:,1:ncat) * Lsub
[3193]	379	ENDIF
[9500]	380
[3625]	381	DO jl=1,ncat
[9819]	382	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1.,'flatn' )
[2874]	383
	384	! GBM conductive flux through ice (CI_6)
	385	! Convert to GBM
[9500]	386	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
[9675]	387	ztmp(:,:) = qcn_ice(:,:,jl)*a_i(:,:,jl)
[3193]	388	ELSE
[9675]	389	ztmp(:,:) = qcn_ice(:,:,jl)
[3193]	390	ENDIF
[9819]	391	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1.,'fcondtopn' )
[2874]	392
	393	! GBM surface heat flux (CI_7)
	394	! Convert to GBM
[9500]	395	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
[9675]	396	ztmp(:,:) = (qml_ice(:,:,jl)+qcn_ice(:,:,jl))*a_i(:,:,jl)
[3193]	397	ELSE
[9675]	398	ztmp(:,:) = (qml_ice(:,:,jl)+qcn_ice(:,:,jl))
[3193]	399	ENDIF
[9819]	400	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. ,'fsurfn_f')
[3193]	401	ENDDO
[2874]	402
[7646]	403	ELSE IF (ksbc == jp_blk) THEN
[2874]	404
[7646]	405	! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
[2874]	406	! x comp and y comp of atmosphere surface wind (CICE expects on T points)
[3193]	407	ztmp(:,:) = wndi_ice(:,:)
[9819]	408	CALL nemo2cice(ztmp,uatm,'T', -1. ,'uatm')
[3193]	409	ztmp(:,:) = wndj_ice(:,:)
[9819]	410	CALL nemo2cice(ztmp,vatm,'T', -1.,'vatm' )
[3193]	411	ztmp(:,:) = SQRT ( wndi_ice(:,:)2 + wndj_ice(:,:)2 )
[9819]	412	CALL nemo2cice(ztmp,wind,'T', 1.,'wind' ) ! Wind speed (m/s)
[3193]	413	ztmp(:,:) = qsr_ice(:,:,1)
[9819]	414	CALL nemo2cice(ztmp,fsw,'T', 1.,'fsw' ) ! Incoming short-wave (W/m^2)
[3193]	415	ztmp(:,:) = qlw_ice(:,:,1)
[9819]	416	CALL nemo2cice(ztmp,flw,'T', 1.,'flw' ) ! Incoming long-wave (W/m^2)
[3193]	417	ztmp(:,:) = tatm_ice(:,:)
[9819]	418	CALL nemo2cice(ztmp,Tair,'T', 1.,'tair' ) ! Air temperature (K)
	419	CALL nemo2cice(ztmp,potT,'T', 1.,'pott' ) ! Potential temp (K)
[2874]	420	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
[3193]	421	ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
	422	! Constant (101000.) atm pressure assumed
[9819]	423	CALL nemo2cice(ztmp,rhoa,'T', 1.,'rhoa' ) ! Air density (kg/m^3)
[3193]	424	ztmp(:,:) = qatm_ice(:,:)
[9819]	425	CALL nemo2cice(ztmp,Qa,'T', 1.,'qa' ) ! Specific humidity (kg/kg)
[3193]	426	ztmp(:,:)=10.0
[9819]	427	CALL nemo2cice(ztmp,zlvl,'T', 1.,'zlvl' ) ! Atmos level height (m)
[2874]	428
	429	! May want to check all values are physically realistic (as in CICE routine
	430	! prepare_forcing)?
	431
	432	! Divide shortwave into spectral bands (as in prepare_forcing)
[3193]	433	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
[9819]	434	CALL nemo2cice(ztmp,swvdr,'T', 1.,'swvdr' )
[3193]	435	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
[9819]	436	CALL nemo2cice(ztmp,swvdf,'T', 1.,'swvdf' )
[3193]	437	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
[9819]	438	CALL nemo2cice(ztmp,swidr,'T', 1.,'swidr' )
[3193]	439	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
[9819]	440	CALL nemo2cice(ztmp,swidf,'T', 1. ,'swidf')
[2874]	441
	442	ENDIF
	443
[9500]	444	#if defined key_asminc
	445	!Ice concentration change (from assimilation)
	446	ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
[9819]	447	Call nemo2cice(ztmp,daice_da,'T', 1. ,'daice_da')
[9500]	448	#endif
	449
[2874]	450	! Snowfall
[4990]	451	! Ensure fsnow is positive (as in CICE routine prepare_forcing)
	452	IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
[3193]	453	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
[9819]	454	CALL nemo2cice(ztmp,fsnow,'T', 1.,'fsnow' )
[2874]	455
	456	! Rainfall
[4990]	457	IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
[3193]	458	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
[9819]	459	CALL nemo2cice(ztmp,frain,'T', 1.,'frain' )
[2874]	460
[9500]	461	! Recalculate freezing temperature and send to CICE
	462	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1))
[9819]	463	CALL nemo2cice(sstfrz,Tf,'T', 1. ,'Tf')
[9500]	464
[2874]	465	! Freezing/melting potential
[3275]	466	! Calculated over NEMO leapfrog timestep (hence 2*dt)
[9500]	467	nfrzmlt(:,:)=rau0rcpe3t_m(:,:)(sstfrz(:,:)-sst_m(:,:))/(2.0dt)
[9819]	468	CALL nemo2cice(nfrzmlt,frzmlt,'T', 1.,'frzmlt' )
[2874]	469
	470	! SST and SSS
	471
[9819]	472	CALL nemo2cice(sst_m,sst,'T', 1.,'sst' )
	473	CALL nemo2cice(sss_m,sss,'T', 1.,'sss' )
[2874]	474
[9500]	475	IF( ksbc == jp_purecpl ) THEN
	476	! Sea ice surface skin temperature
	477	DO jl=1,ncat
[9819]	478	CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1., 'trcrn')
[9500]	479	ENDDO
	480	ENDIF
	481
[2874]	482	! x comp and y comp of surface ocean current
	483	! U point to F point
[3193]	484	DO jj=1,jpjm1
	485	DO ji=1,jpi
	486	ztmp(ji,jj)=0.5(ssu_m(ji,jj)+ssu_m(ji,jj+1))fmask(ji,jj,1)
	487	ENDDO
	488	ENDDO
[9819]	489	CALL nemo2cice(ztmp,uocn,'F', -1., 'uocn' )
[2874]	490
	491	! V point to F point
[3193]	492	DO jj=1,jpj
	493	DO ji=1,jpim1
	494	ztmp(ji,jj)=0.5(ssv_m(ji,jj)+ssv_m(ji+1,jj))fmask(ji,jj,1)
	495	ENDDO
	496	ENDDO
[9819]	497	CALL nemo2cice(ztmp,vocn,'F', -1. ,'vocn')
[2874]	498
[9019]	499	IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==!
[3625]	500	!
	501	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	502	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
	503	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
	504	!
	505	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	506	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
	507	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
	508	!
	509	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
	510	!
	511	!
	512	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	513	zpice(:,:) = ssh_m(:,:)
	514	ENDIF
	515
[3189]	516	! x comp and y comp of sea surface slope (on F points)
	517	! T point to F point
[5836]	518	DO jj = 1, jpjm1
	519	DO ji = 1, jpim1
	520	ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj )) * r1_e1u(ji,jj ) &
	521	& + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1) ) * fmask(ji,jj,1)
	522	END DO
	523	END DO
[9819]	524	CALL nemo2cice( ztmp,ss_tltx,'F', -1. , 'ss_tltx')
[3189]	525
	526	! T point to F point
[5836]	527	DO jj = 1, jpjm1
	528	DO ji = 1, jpim1
	529	ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj)) * r1_e2v(ji ,jj) &
	530	& + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj) ) * fmask(ji,jj,1)
	531	END DO
	532	END DO
[9819]	533	CALL nemo2cice(ztmp,ss_tlty,'F', -1. ,'ss_tlty')
[3193]	534	!
[2874]	535	END SUBROUTINE cice_sbc_in
	536
[3152]	537
[5836]	538	SUBROUTINE cice_sbc_out( kt, ksbc )
[2874]	539	!!---------------------------------------------------------------------
	540	!! * ROUTINE cice_sbc_out *
[3040]	541	!! ** Purpose: Get fields from CICE and set surface fields for NEMO
[3152]	542	!!---------------------------------------------------------------------
[2874]	543	INTEGER, INTENT( in ) :: kt ! ocean time step
[4990]	544	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
[3152]	545
[3625]	546	INTEGER :: ji, jj, jl ! dummy loop indices
[9125]	547	REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
[2874]	548	!!---------------------------------------------------------------------
[3193]	549	!
[3152]	550	IF( kt == nit000 ) THEN
[2874]	551	IF(lwp) WRITE(numout,*)'cice_sbc_out'
[3152]	552	ENDIF
	553
[2874]	554	! x comp of ocean-ice stress
[3625]	555	CALL cice2nemo(strocnx,ztmp1,'F', -1. )
[3193]	556	ss_iou(:,:)=0.0
[2874]	557	! F point to U point
[3193]	558	DO jj=2,jpjm1
	559	DO ji=2,jpim1
[3625]	560	ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
[3193]	561	ENDDO
	562	ENDDO
	563	CALL lbc_lnk( ss_iou , 'U', -1. )
[2874]	564
	565	! y comp of ocean-ice stress
[3625]	566	CALL cice2nemo(strocny,ztmp1,'F', -1. )
[3193]	567	ss_iov(:,:)=0.0
[2874]	568	! F point to V point
	569
[3193]	570	DO jj=1,jpjm1
	571	DO ji=2,jpim1
[3625]	572	ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
[3193]	573	ENDDO
	574	ENDDO
	575	CALL lbc_lnk( ss_iov , 'V', -1. )
[2874]	576
	577	! x and y comps of surface stress
	578	! Combine wind stress and ocean-ice stress
	579	! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
[5133]	580	! strocnx and strocny already weighted by ice fraction in CICE so not done here
[2874]	581
[3193]	582	utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
	583	vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
[5133]	584
	585	! Also need ice/ocean stress on T points so that taum can be updated
	586	! This interpolation is already done in CICE so best to use those values
	587	CALL cice2nemo(strocnxT,ztmp1,'T',-1.)
	588	CALL cice2nemo(strocnyT,ztmp2,'T',-1.)
	589
	590	! Update taum with modulus of ice-ocean stress
	591	! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
[5836]	592	taum(:,:)=(1.0-fr_i(:,:))taum(:,:)+fr_i(:,:)SQRT(ztmp1ztmp1 + ztmp2ztmp2)
[2874]	593
	594	! Freshwater fluxes
	595
[4990]	596	IF (ksbc == jp_flx) THEN
[2874]	597	! Note that emp from the forcing files is evap(1-aice)-(tprecip-aicesprecip)
	598	! What we want here is evap(1-aice)-tprecip(1-aice) hence manipulation below
	599	! Not ideal since aice won't be the same as in the atmosphere.
	600	! Better to use evap and tprecip? (but for now don't read in evap in this case)
[3193]	601	emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
[7646]	602	ELSE IF (ksbc == jp_blk) THEN
[3193]	603	emp(:,:) = (1.0-fr_i(:,:))*emp(:,:)
[5407]	604	ELSE IF (ksbc == jp_purecpl) THEN
[3625]	605	! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
	606	! This is currently as required with the coupling fields from the UM atmosphere
[3193]	607	emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
	608	ENDIF
[2874]	609
[4990]	610	#if defined key_cice4
[3625]	611	CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
	612	CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
[4990]	613	#else
	614	CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
	615	CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
	616	#endif
[2874]	617
[3625]	618	! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
	619	! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
	620	! which has to be compensated for by the ocean salinity potentially going negative
	621	! This check breaks conservation but seems reasonable until we have prognostic ice salinity
	622	! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
	623	WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
	624	sfx(:,:)=ztmp2(:,:)*1000.0
	625	emp(:,:)=emp(:,:)-ztmp1(:,:)
[4990]	626	fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
	627
[9098]	628	CALL lbc_lnk_multi( emp , 'T', 1., sfx , 'T', 1. )
[2874]	629
	630	! Solar penetrative radiation and non solar surface heat flux
	631
	632	! Scale qsr and qns according to ice fraction (bulk formulae only)
	633
[7646]	634	IF (ksbc == jp_blk) THEN
[3193]	635	qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
	636	qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
	637	ENDIF
[2874]	638	! Take into account snow melting except for fully coupled when already in qns_tot
[5407]	639	IF (ksbc == jp_purecpl) THEN
[3193]	640	qsr(:,:)= qsr_tot(:,:)
	641	qns(:,:)= qns_tot(:,:)
	642	ELSE
	643	qns(:,:)= qns(:,:)-sprecip(:,:)Lfresh(1.0-fr_i(:,:))
	644	ENDIF
[2874]	645
	646	! Now add in ice / snow related terms
	647	! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
[4990]	648	#if defined key_cice4
[3625]	649	CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
[4990]	650	#else
	651	CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
	652	#endif
[3625]	653	qsr(:,:)=qsr(:,:)+ztmp1(:,:)
[3193]	654	CALL lbc_lnk( qsr , 'T', 1. )
[2874]	655
[3193]	656	DO jj=1,jpj
	657	DO ji=1,jpi
[2874]	658	nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
[3193]	659	ENDDO
	660	ENDDO
[2874]	661
[4990]	662	#if defined key_cice4
[3625]	663	CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
[4990]	664	#else
	665	CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
	666	#endif
[3625]	667	qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
[2874]	668
[3193]	669	CALL lbc_lnk( qns , 'T', 1. )
[2874]	670
	671	! Prepare for the following CICE time-step
	672
[3193]	673	CALL cice2nemo(aice,fr_i,'T', 1. )
[5407]	674	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
[3625]	675	DO jl=1,ncat
	676	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	677	ENDDO
	678	ENDIF
[2874]	679
	680	! T point to U point
	681	! T point to V point
[3193]	682	DO jj=1,jpjm1
	683	DO ji=1,jpim1
	684	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	685	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	686	ENDDO
	687	ENDDO
[2874]	688
[9098]	689	CALL lbc_lnk_multi( fr_iu , 'U', 1., fr_iv , 'V', 1. )
[2874]	690
[9019]	691	! set the snow+ice mass
	692	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	693	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	694	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	695	snwice_mass_b(:,:) = snwice_mass(:,:)
	696	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
[3193]	697	!
[2874]	698	END SUBROUTINE cice_sbc_out
	699
[3152]	700
[2874]	701	SUBROUTINE cice_sbc_hadgam( kt )
	702	!!---------------------------------------------------------------------
	703	!! * ROUTINE cice_sbc_hadgam *
[3040]	704	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
[2874]	705	!!
	706	!!
[9124]	707	!!---------------------------------------------------------------------
[2874]	708	INTEGER, INTENT( in ) :: kt ! ocean time step
[9124]	709	!!
[3625]	710	INTEGER :: jl ! dummy loop index
[3193]	711	INTEGER :: ierror
[9124]	712	!!---------------------------------------------------------------------
[9500]	713	#if defined key_asminc
	714	! Import fresh water and salt flux due to seaice da
	715	CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
	716	CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
	717	#endif
	718
[9498]	719	!
[2874]	720	! ! =========================== !
	721	! ! Prepare Coupling fields !
	722	! ! =========================== !
[9019]	723	!
	724	! x and y comp of ice velocity
	725	!
[3193]	726	CALL cice2nemo(uvel,u_ice,'F', -1. )
	727	CALL cice2nemo(vvel,v_ice,'F', -1. )
	728	!
[9019]	729	! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
	730	!
	731	! Snow and ice thicknesses (CO_2 and CO_3)
	732	!
	733	DO jl = 1, ncat
	734	CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
	735	CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
	736	END DO
[9500]	737
	738	#if ! defined key_cice4
	739	! Meltpond fraction and depth
	740	DO jl = 1,ncat
[9675]	741	CALL cice2nemo(apeffn(:,:,jl,:),a_ip(:,:,jl),'T', 1. )
	742	CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),v_ip(:,:,jl),'T', 1. )
[9500]	743	ENDDO
[9695]	744
	745	v_ip(:,:,:) = v_ip(:,:,:) * a_i(:,:,:)
[9500]	746	#endif
	747
	748
	749	! If using multilayers thermodynamics in CICE then get top layer temperature
	750	! and effective conductivity
	751	!! When using NEMO with CICE, this change requires use of
	752	!! one of the following two CICE branches:
	753	!! - at CICE5.0, hadax/r1015_GSI8_with_GSI7
	754	!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
	755	IF (heat_capacity) THEN
	756	DO jl = 1,ncat
[9675]	757	CALL cice2nemo(Tn_top(:,:,jl,:),t1_ice(:,:,jl),'T', 1. )
	758	CALL cice2nemo(keffn_top(:,:,jl,:),cnd_ice(:,:,jl),'T', 1. )
[9500]	759	ENDDO
	760	! Convert surface temperature to Kelvin
[9675]	761	t1_ice(:,:,:)=t1_ice(:,:,:)+rt0
[9500]	762	ELSE
[9675]	763	t1_ice(:,:,:) = 0.0
	764	cnd_ice(:,:,:) = 0.0
[9500]	765	ENDIF
	766
[9019]	767	!
[2874]	768	END SUBROUTINE cice_sbc_hadgam
	769
	770
	771	SUBROUTINE cice_sbc_final
	772	!!---------------------------------------------------------------------
	773	!! * ROUTINE cice_sbc_final *
	774	!! ** Purpose: Finalize CICE
	775	!!---------------------------------------------------------------------
[9124]	776	!
[2874]	777	IF(lwp) WRITE(numout,*)'cice_sbc_final'
[9124]	778	!
[3193]	779	CALL CICE_Finalize
[9124]	780	!
[2874]	781	END SUBROUTINE cice_sbc_final
	782
[9124]	783
[2874]	784	SUBROUTINE cice_sbc_force (kt)
	785	!!---------------------------------------------------------------------
	786	!! * ROUTINE cice_sbc_force *
	787	!! ** Purpose : Provide CICE forcing from files
	788	!!
	789	!!---------------------------------------------------------------------
	790	!! ** Method : READ monthly flux file in NetCDF files
	791	!!
	792	!! snowfall
	793	!! rainfall
	794	!! sublimation rate
	795	!! topmelt (category)
	796	!! botmelt (category)
	797	!!
	798	!! History :
	799	!!----------------------------------------------------------------------
	800	USE iom
[9124]	801	!!
	802	INTEGER, INTENT( in ) :: kt ! ocean time step
	803	!!
[2874]	804	INTEGER :: ierror ! return error code
	805	INTEGER :: ifpr ! dummy loop index
	806	!!
	807	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
	808	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
	809	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
	810	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
	811	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
	812	!!
	813	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
	814	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
	815	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
[4230]	816	INTEGER :: ios
[2874]	817	!!---------------------------------------------------------------------
	818
	819	! ! ====================== !
	820	IF( kt == nit000 ) THEN ! First call kt=nit000 !
	821	! ! ====================== !
[4990]	822	! namsbc_cice is not yet in the reference namelist
	823	! set file information (default values)
	824	cn_dir = './' ! directory in which the model is executed
	825
	826	! (NB: frequency positive => hours, negative => months)
	827	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask
	828	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file
	829	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	830	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	831	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' )
	832	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	833	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	834	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	835	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	836	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	837	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	838	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	839	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	840	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	841	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	842
[4230]	843	REWIND( numnam_ref ) ! Namelist namsbc_cice in reference namelist :
	844	READ ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
[9168]	845	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
[2874]	846
[4230]	847	REWIND( numnam_cfg ) ! Namelist namsbc_cice in configuration namelist : Parameters of the run
	848	READ ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
[9168]	849	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
[4624]	850	IF(lwm) WRITE ( numond, namsbc_cice )
[2874]	851
	852	! store namelist information in an array
	853	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
	854	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
	855	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
	856	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
	857	slf_i(jp_bot5) = sn_bot5
	858
	859	! set sf structure
	860	ALLOCATE( sf(jpfld), STAT=ierror )
	861	IF( ierror > 0 ) THEN
	862	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
	863	ENDIF
	864
	865	DO ifpr= 1, jpfld
	866	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
	867	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
	868	END DO
	869
	870	! fill sf with slf_i and control print
	871	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
	872	!
	873	ENDIF
	874
	875	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
	876	! ! input fields at the current time-step
	877
	878	! set the fluxes from read fields
	879	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
	880	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
[3040]	881	! May be better to do this conversion somewhere else
[9675]	882	zevap_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
	883	qml_ice(:,:,1) = sf(jp_top1)%fnow(:,:,1)
	884	qml_ice(:,:,2) = sf(jp_top2)%fnow(:,:,1)
	885	qml_ice(:,:,3) = sf(jp_top3)%fnow(:,:,1)
	886	qml_ice(:,:,4) = sf(jp_top4)%fnow(:,:,1)
	887	qml_ice(:,:,5) = sf(jp_top5)%fnow(:,:,1)
	888	qcn_ice(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
	889	qcn_ice(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
	890	qcn_ice(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
	891	qcn_ice(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
	892	qcn_ice(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
[2874]	893
	894	! control print (if less than 100 time-step asked)
	895	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
	896	WRITE(numout,*)
	897	WRITE(numout,*) ' read forcing fluxes for CICE OK'
	898	CALL FLUSH(numout)
	899	ENDIF
	900
	901	END SUBROUTINE cice_sbc_force
	902
[9819]	903	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn, varname)
[2874]	904	!!---------------------------------------------------------------------
	905	!! * ROUTINE nemo2cice *
	906	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
	907	#if defined key_nemocice_decomp
	908	!!
	909	!! NEMO and CICE PE sub domains are identical, hence
	910	!! there is no need to gather or scatter data from
	911	!! one PE configuration to another.
	912	#else
	913	!! Automatically gather/scatter between
	914	!! different processors and blocks
	915	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
	916	!! B. Gather pn into global array (png)
	917	!! C. Map png into CICE global array (pcg)
	918	!! D. Scatter pcg to CICE blocks (pc) + update haloes
	919	#endif
	920	!!---------------------------------------------------------------------
[3193]	921	CHARACTER(len=1), INTENT( in ) :: &
	922	cd_type ! nature of pn grid-point
	923	! ! = T or F gridpoints
	924	REAL(wp), INTENT( in ) :: &
	925	psgn ! control of the sign change
	926	! ! =-1 , the sign is modified following the type of b.c. used
	927	! ! = 1 , no sign change
	928	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	929	#if !defined key_nemocice_decomp
[3625]	930	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
[3193]	931	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	932	#endif
[3193]	933	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
[9819]	934	CHARACTER(len=10), INTENT( in ) :: varname
[3193]	935	INTEGER (int_kind) :: &
	936	field_type, &! id for type of field (scalar, vector, angle)
	937	grid_loc ! id for location on horizontal grid
[2874]	938	! (center, NEcorner, Nface, Eface)
	939
[3193]	940	INTEGER :: ji, jj, jn ! dummy loop indices
[5836]	941	!!---------------------------------------------------------------------
[2874]	942
[3193]	943	! A. Ensure all haloes are filled in NEMO field (pn)
[2874]	944
[3193]	945	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	946
	947	#if defined key_nemocice_decomp
	948
[3193]	949	! Copy local domain data from NEMO to CICE field
	950	pc(:,:,1)=0.0
[9631]	951	DO jj=2,jpj-1
	952	DO ji=2,jpi-1
[3625]	953	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
[3193]	954	ENDDO
	955	ENDDO
[2874]	956
	957	#else
	958
[3193]	959	! B. Gather pn into global array (png)
[2874]	960
[3193]	961	IF ( jpnij > 1) THEN
	962	CALL mppsync
	963	CALL mppgather (pn,0,png)
	964	CALL mppsync
	965	ELSE
	966	png(:,:,1)=pn(:,:)
	967	ENDIF
[2874]	968
[3193]	969	! C. Map png into CICE global array (pcg)
[2874]	970
	971	! Need to make sure this is robust to changes in NEMO halo rows....
	972	! (may be OK but not 100% sure)
	973
[3193]	974	IF (nproc==0) THEN
[2874]	975	! pcg(:,:)=0.0
[3193]	976	DO jn=1,jpnij
[3625]	977	DO jj=nldjt(jn),nlejt(jn)
	978	DO ji=nldit(jn),nleit(jn)
	979	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
[3193]	980	ENDDO
	981	ENDDO
	982	ENDDO
[3625]	983	DO jj=1,ny_global
	984	DO ji=1,nx_global
	985	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
	986	ENDDO
	987	ENDDO
[3193]	988	ENDIF
[2874]	989
	990	#endif
	991
[3193]	992	SELECT CASE ( cd_type )
	993	CASE ( 'T' )
	994	grid_loc=field_loc_center
	995	CASE ( 'F' )
	996	grid_loc=field_loc_NEcorner
	997	END SELECT
[2874]	998
[3193]	999	SELECT CASE ( NINT(psgn) )
	1000	CASE ( -1 )
	1001	field_type=field_type_vector
	1002	CASE ( 1 )
	1003	field_type=field_type_scalar
	1004	END SELECT
[2874]	1005
	1006	#if defined key_nemocice_decomp
[3193]	1007	! Ensure CICE halos are up to date
	1008	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1009	#else
[3193]	1010	! D. Scatter pcg to CICE blocks (pc) + update halos
	1011	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
[2874]	1012	#endif
	1013
[9819]	1014	IF ( ln_ctl ) THEN
	1015	WRITE(numout,*)'nemo2cice: ',varname,' min,max = ',MINVAL(pc(:,:,:)), MAXVAL(pc(:,:,:))
	1016	CALL flush(numout)
	1017	ENDIF
[2874]	1018	END SUBROUTINE nemo2cice
	1019
	1020	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
	1021	!!---------------------------------------------------------------------
	1022	!! * ROUTINE cice2nemo *
	1023	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
	1024	#if defined key_nemocice_decomp
	1025	!!
	1026	!! NEMO and CICE PE sub domains are identical, hence
	1027	!! there is no need to gather or scatter data from
	1028	!! one PE configuration to another.
	1029	#else
	1030	!! Automatically deal with scatter/gather between
	1031	!! different processors and blocks
	1032	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
	1033	!! B. Map pcg into NEMO global array (png)
	1034	!! C. Scatter png into NEMO field (pn) for each processor
	1035	!! D. Ensure all haloes are filled in pn
	1036	#endif
	1037	!!---------------------------------------------------------------------
	1038
[3193]	1039	CHARACTER(len=1), INTENT( in ) :: &
	1040	cd_type ! nature of pn grid-point
	1041	! ! = T or F gridpoints
	1042	REAL(wp), INTENT( in ) :: &
	1043	psgn ! control of the sign change
	1044	! ! =-1 , the sign is modified following the type of b.c. used
	1045	! ! = 1 , no sign change
	1046	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	1047
	1048	#if defined key_nemocice_decomp
[3193]	1049	INTEGER (int_kind) :: &
	1050	field_type, & ! id for type of field (scalar, vector, angle)
	1051	grid_loc ! id for location on horizontal grid
	1052	! (center, NEcorner, Nface, Eface)
[2874]	1053	#else
[3193]	1054	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	1055	#endif
	1056
[3193]	1057	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
[2874]	1058
[3193]	1059	INTEGER :: ji, jj, jn ! dummy loop indices
[2874]	1060
	1061
	1062	#if defined key_nemocice_decomp
	1063
[3193]	1064	SELECT CASE ( cd_type )
	1065	CASE ( 'T' )
	1066	grid_loc=field_loc_center
	1067	CASE ( 'F' )
	1068	grid_loc=field_loc_NEcorner
	1069	END SELECT
[2874]	1070
[3193]	1071	SELECT CASE ( NINT(psgn) )
	1072	CASE ( -1 )
	1073	field_type=field_type_vector
	1074	CASE ( 1 )
	1075	field_type=field_type_scalar
	1076	END SELECT
[2874]	1077
[3193]	1078	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1079
	1080
[3193]	1081	pn(:,:)=0.0
	1082	DO jj=1,jpjm1
	1083	DO ji=1,jpim1
[3625]	1084	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
[3193]	1085	ENDDO
	1086	ENDDO
[2874]	1087
	1088	#else
	1089
[3193]	1090	! A. Gather CICE blocks (pc) into global array (pcg)
[2874]	1091
[3193]	1092	CALL gather_global(pcg, pc, 0, distrb_info)
[2874]	1093
	1094	! B. Map pcg into NEMO global array (png)
	1095
	1096	! Need to make sure this is robust to changes in NEMO halo rows....
	1097	! (may be OK but not spent much time thinking about it)
[3625]	1098	! Note that non-existent pcg elements may be used below, but
	1099	! the lbclnk call on pn will replace these with sensible values
[2874]	1100
[3193]	1101	IF (nproc==0) THEN
	1102	png(:,:,:)=0.0
	1103	DO jn=1,jpnij
[3625]	1104	DO jj=nldjt(jn),nlejt(jn)
	1105	DO ji=nldit(jn),nleit(jn)
	1106	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
[3193]	1107	ENDDO
	1108	ENDDO
	1109	ENDDO
	1110	ENDIF
[2874]	1111
[3193]	1112	! C. Scatter png into NEMO field (pn) for each processor
[2874]	1113
[3193]	1114	IF ( jpnij > 1) THEN
	1115	CALL mppsync
	1116	CALL mppscatter (png,0,pn)
	1117	CALL mppsync
	1118	ELSE
	1119	pn(:,:)=png(:,:,1)
	1120	ENDIF
[2874]	1121
	1122	#endif
	1123
[3193]	1124	! D. Ensure all haloes are filled in pn
[2874]	1125
[3193]	1126	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	1127
	1128	END SUBROUTINE cice2nemo
	1129
	1130	#else
	1131	!!----------------------------------------------------------------------
	1132	!! Default option Dummy module NO CICE sea-ice model
	1133	!!----------------------------------------------------------------------
	1134	CONTAINS
	1135
[4990]	1136	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
[2874]	1137	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
	1138	END SUBROUTINE sbc_ice_cice
	1139
[4990]	1140	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
[2874]	1141	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
	1142	END SUBROUTINE cice_sbc_init
	1143
	1144	SUBROUTINE cice_sbc_final ! Dummy routine
	1145	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
	1146	END SUBROUTINE cice_sbc_final
	1147
	1148	#endif
	1149
	1150	!!======================================================================
	1151	END MODULE sbcice_cice

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

Download in other formats: