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.

sbcice_cice.F90 in branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

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

Last change on this file since 9500 was 9500, checked in by davestorkey, 6 years ago
branches/UKMO/dev_merge_2017_CICE_interface : recommit science changes.
File size: 46.2 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
	211	CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)
	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))
	218	CALL nemo2cice(sstfrz,Tf, 'T', 1. )
	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
	342	CALL nemo2cice(ztmp,strax,'F', -1. )
[2874]	343
	344	! y comp of wind stress (CI_2)
	345	! V point to F point
[3193]	346	DO jj=1,jpj
	347	DO ji=1,jpim1
	348	ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
	349	+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
	350	ENDDO
	351	ENDDO
	352	CALL nemo2cice(ztmp,stray,'F', -1. )
[2874]	353
[9500]	354
	355	! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in
	356	! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby
	357	! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing
	358	! gridbox mean fluxes in the UM by future ice concentration obtained through
	359	! OASIS. This allows for a much more realistic apportionment of energy through
	360	! the ice - and conserves energy.
	361	! Therefore the fluxes are now divided by ice concentration in the coupled
	362	! formulation (jp_purecpl) as well as for jp_flx. This NEMO branch should only
	363	! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at
	364	! which point the GSI8 UM changes were committed.
	365
[2874]	366	! Surface downward latent heat flux (CI_5)
[4990]	367	IF (ksbc == jp_flx) THEN
[3625]	368	DO jl=1,ncat
	369	ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
[3193]	370	ENDDO
[9500]	371	ELSE IF (ksbc == jp_purecpl) THEN
	372	DO jl=1,ncat
	373	ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl)
[3193]	374	ENDDO
[9500]	375	ELSE
	376	!In coupled mode - qla_ice calculated in sbc_cpl for each category
	377	ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat)
[3193]	378	ENDIF
[9500]	379
[3625]	380	DO jl=1,ncat
	381	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
[2874]	382
	383	! GBM conductive flux through ice (CI_6)
	384	! Convert to GBM
[9500]	385	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
[3625]	386	ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
[3193]	387	ELSE
[3625]	388	ztmp(:,:) = botmelt(:,:,jl)
[3193]	389	ENDIF
[3625]	390	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
[2874]	391
	392	! GBM surface heat flux (CI_7)
	393	! Convert to GBM
[9500]	394	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
[3625]	395	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)
[3193]	396	ELSE
[3625]	397	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
[3193]	398	ENDIF
[3625]	399	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
[3193]	400	ENDDO
[2874]	401
[7646]	402	ELSE IF (ksbc == jp_blk) THEN
[2874]	403
[7646]	404	! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
[2874]	405	! x comp and y comp of atmosphere surface wind (CICE expects on T points)
[3193]	406	ztmp(:,:) = wndi_ice(:,:)
	407	CALL nemo2cice(ztmp,uatm,'T', -1. )
	408	ztmp(:,:) = wndj_ice(:,:)
	409	CALL nemo2cice(ztmp,vatm,'T', -1. )
	410	ztmp(:,:) = SQRT ( wndi_ice(:,:)2 + wndj_ice(:,:)2 )
	411	CALL nemo2cice(ztmp,wind,'T', 1. ) ! Wind speed (m/s)
	412	ztmp(:,:) = qsr_ice(:,:,1)
	413	CALL nemo2cice(ztmp,fsw,'T', 1. ) ! Incoming short-wave (W/m^2)
	414	ztmp(:,:) = qlw_ice(:,:,1)
	415	CALL nemo2cice(ztmp,flw,'T', 1. ) ! Incoming long-wave (W/m^2)
	416	ztmp(:,:) = tatm_ice(:,:)
	417	CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K)
	418	CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K)
[2874]	419	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
[3193]	420	ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
	421	! Constant (101000.) atm pressure assumed
	422	CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3)
	423	ztmp(:,:) = qatm_ice(:,:)
	424	CALL nemo2cice(ztmp,Qa,'T', 1. ) ! Specific humidity (kg/kg)
	425	ztmp(:,:)=10.0
	426	CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m)
[2874]	427
	428	! May want to check all values are physically realistic (as in CICE routine
	429	! prepare_forcing)?
	430
	431	! Divide shortwave into spectral bands (as in prepare_forcing)
[3193]	432	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
[2874]	433	CALL nemo2cice(ztmp,swvdr,'T', 1. )
[3193]	434	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
[2874]	435	CALL nemo2cice(ztmp,swvdf,'T', 1. )
[3193]	436	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
[2874]	437	CALL nemo2cice(ztmp,swidr,'T', 1. )
[3193]	438	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
[2874]	439	CALL nemo2cice(ztmp,swidf,'T', 1. )
	440
	441	ENDIF
	442
[9500]	443	#if defined key_asminc
	444	!Ice concentration change (from assimilation)
	445	ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
	446	Call nemo2cice(ztmp,daice_da,'T', 1. )
	447	#endif
	448
[2874]	449	! Snowfall
[4990]	450	! Ensure fsnow is positive (as in CICE routine prepare_forcing)
	451	IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
[3193]	452	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
	453	CALL nemo2cice(ztmp,fsnow,'T', 1. )
[2874]	454
	455	! Rainfall
[4990]	456	IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
[3193]	457	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
	458	CALL nemo2cice(ztmp,frain,'T', 1. )
[2874]	459
[9500]	460	! Recalculate freezing temperature and send to CICE
	461	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1))
	462	CALL nemo2cice(sstfrz,Tf,'T', 1. )
	463
[2874]	464	! Freezing/melting potential
[3275]	465	! Calculated over NEMO leapfrog timestep (hence 2*dt)
[9500]	466	nfrzmlt(:,:)=rau0rcpe3t_m(:,:)(sstfrz(:,:)-sst_m(:,:))/(2.0dt)
	467	CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. )
[2874]	468
	469	! SST and SSS
	470
[3193]	471	CALL nemo2cice(sst_m,sst,'T', 1. )
	472	CALL nemo2cice(sss_m,sss,'T', 1. )
[2874]	473
[9500]	474	IF( ksbc == jp_purecpl ) THEN
	475	! Sea ice surface skin temperature
	476	DO jl=1,ncat
	477	CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.)
	478	ENDDO
	479	ENDIF
	480
[2874]	481	! x comp and y comp of surface ocean current
	482	! U point to F point
[3193]	483	DO jj=1,jpjm1
	484	DO ji=1,jpi
	485	ztmp(ji,jj)=0.5(ssu_m(ji,jj)+ssu_m(ji,jj+1))fmask(ji,jj,1)
	486	ENDDO
	487	ENDDO
	488	CALL nemo2cice(ztmp,uocn,'F', -1. )
[2874]	489
	490	! V point to F point
[3193]	491	DO jj=1,jpj
	492	DO ji=1,jpim1
	493	ztmp(ji,jj)=0.5(ssv_m(ji,jj)+ssv_m(ji+1,jj))fmask(ji,jj,1)
	494	ENDDO
	495	ENDDO
	496	CALL nemo2cice(ztmp,vocn,'F', -1. )
[2874]	497
[9019]	498	IF( ln_ice_embd ) THEN !== embedded sea ice: compute representative ice top surface ==!
[3625]	499	!
	500	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	501	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
	502	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
	503	!
	504	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	505	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
	506	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
	507	!
	508	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
	509	!
	510	!
	511	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	512	zpice(:,:) = ssh_m(:,:)
	513	ENDIF
	514
[3189]	515	! x comp and y comp of sea surface slope (on F points)
	516	! T point to F point
[5836]	517	DO jj = 1, jpjm1
	518	DO ji = 1, jpim1
	519	ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj )) * r1_e1u(ji,jj ) &
	520	& + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1) ) * fmask(ji,jj,1)
	521	END DO
	522	END DO
	523	CALL nemo2cice( ztmp,ss_tltx,'F', -1. )
[3189]	524
	525	! T point to F point
[5836]	526	DO jj = 1, jpjm1
	527	DO ji = 1, jpim1
	528	ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj)) * r1_e2v(ji ,jj) &
	529	& + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj) ) * fmask(ji,jj,1)
	530	END DO
	531	END DO
[3193]	532	CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
	533	!
[2874]	534	END SUBROUTINE cice_sbc_in
	535
[3152]	536
[5836]	537	SUBROUTINE cice_sbc_out( kt, ksbc )
[2874]	538	!!---------------------------------------------------------------------
	539	!! * ROUTINE cice_sbc_out *
[3040]	540	!! ** Purpose: Get fields from CICE and set surface fields for NEMO
[3152]	541	!!---------------------------------------------------------------------
[2874]	542	INTEGER, INTENT( in ) :: kt ! ocean time step
[4990]	543	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
[3152]	544
[3625]	545	INTEGER :: ji, jj, jl ! dummy loop indices
[9125]	546	REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
[2874]	547	!!---------------------------------------------------------------------
[3193]	548	!
[3152]	549	IF( kt == nit000 ) THEN
[2874]	550	IF(lwp) WRITE(numout,*)'cice_sbc_out'
[3152]	551	ENDIF
	552
[2874]	553	! x comp of ocean-ice stress
[3625]	554	CALL cice2nemo(strocnx,ztmp1,'F', -1. )
[3193]	555	ss_iou(:,:)=0.0
[2874]	556	! F point to U point
[3193]	557	DO jj=2,jpjm1
	558	DO ji=2,jpim1
[3625]	559	ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
[3193]	560	ENDDO
	561	ENDDO
	562	CALL lbc_lnk( ss_iou , 'U', -1. )
[2874]	563
	564	! y comp of ocean-ice stress
[3625]	565	CALL cice2nemo(strocny,ztmp1,'F', -1. )
[3193]	566	ss_iov(:,:)=0.0
[2874]	567	! F point to V point
	568
[3193]	569	DO jj=1,jpjm1
	570	DO ji=2,jpim1
[3625]	571	ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
[3193]	572	ENDDO
	573	ENDDO
	574	CALL lbc_lnk( ss_iov , 'V', -1. )
[2874]	575
	576	! x and y comps of surface stress
	577	! Combine wind stress and ocean-ice stress
	578	! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
[5133]	579	! strocnx and strocny already weighted by ice fraction in CICE so not done here
[2874]	580
[3193]	581	utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
	582	vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
[5133]	583
	584	! Also need ice/ocean stress on T points so that taum can be updated
	585	! This interpolation is already done in CICE so best to use those values
	586	CALL cice2nemo(strocnxT,ztmp1,'T',-1.)
	587	CALL cice2nemo(strocnyT,ztmp2,'T',-1.)
	588
	589	! Update taum with modulus of ice-ocean stress
	590	! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
[5836]	591	taum(:,:)=(1.0-fr_i(:,:))taum(:,:)+fr_i(:,:)SQRT(ztmp1ztmp1 + ztmp2ztmp2)
[2874]	592
	593	! Freshwater fluxes
	594
[4990]	595	IF (ksbc == jp_flx) THEN
[2874]	596	! Note that emp from the forcing files is evap(1-aice)-(tprecip-aicesprecip)
	597	! What we want here is evap(1-aice)-tprecip(1-aice) hence manipulation below
	598	! Not ideal since aice won't be the same as in the atmosphere.
	599	! Better to use evap and tprecip? (but for now don't read in evap in this case)
[3193]	600	emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
[7646]	601	ELSE IF (ksbc == jp_blk) THEN
[3193]	602	emp(:,:) = (1.0-fr_i(:,:))*emp(:,:)
[5407]	603	ELSE IF (ksbc == jp_purecpl) THEN
[3625]	604	! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
	605	! This is currently as required with the coupling fields from the UM atmosphere
[3193]	606	emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
	607	ENDIF
[2874]	608
[4990]	609	#if defined key_cice4
[3625]	610	CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
	611	CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
[4990]	612	#else
	613	CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
	614	CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
	615	#endif
[2874]	616
[3625]	617	! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
	618	! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
	619	! which has to be compensated for by the ocean salinity potentially going negative
	620	! This check breaks conservation but seems reasonable until we have prognostic ice salinity
	621	! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
	622	WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
	623	sfx(:,:)=ztmp2(:,:)*1000.0
	624	emp(:,:)=emp(:,:)-ztmp1(:,:)
[4990]	625	fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
	626
[9098]	627	CALL lbc_lnk_multi( emp , 'T', 1., sfx , 'T', 1. )
[2874]	628
	629	! Solar penetrative radiation and non solar surface heat flux
	630
	631	! Scale qsr and qns according to ice fraction (bulk formulae only)
	632
[7646]	633	IF (ksbc == jp_blk) THEN
[3193]	634	qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
	635	qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
	636	ENDIF
[2874]	637	! Take into account snow melting except for fully coupled when already in qns_tot
[5407]	638	IF (ksbc == jp_purecpl) THEN
[3193]	639	qsr(:,:)= qsr_tot(:,:)
	640	qns(:,:)= qns_tot(:,:)
	641	ELSE
	642	qns(:,:)= qns(:,:)-sprecip(:,:)Lfresh(1.0-fr_i(:,:))
	643	ENDIF
[2874]	644
	645	! Now add in ice / snow related terms
	646	! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
[4990]	647	#if defined key_cice4
[3625]	648	CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
[4990]	649	#else
	650	CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
	651	#endif
[3625]	652	qsr(:,:)=qsr(:,:)+ztmp1(:,:)
[3193]	653	CALL lbc_lnk( qsr , 'T', 1. )
[2874]	654
[3193]	655	DO jj=1,jpj
	656	DO ji=1,jpi
[2874]	657	nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
[3193]	658	ENDDO
	659	ENDDO
[2874]	660
[4990]	661	#if defined key_cice4
[3625]	662	CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
[4990]	663	#else
	664	CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
	665	#endif
[3625]	666	qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
[2874]	667
[3193]	668	CALL lbc_lnk( qns , 'T', 1. )
[2874]	669
	670	! Prepare for the following CICE time-step
	671
[3193]	672	CALL cice2nemo(aice,fr_i,'T', 1. )
[5407]	673	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
[3625]	674	DO jl=1,ncat
	675	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	676	ENDDO
	677	ENDIF
[2874]	678
	679	! T point to U point
	680	! T point to V point
[3193]	681	DO jj=1,jpjm1
	682	DO ji=1,jpim1
	683	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	684	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	685	ENDDO
	686	ENDDO
[2874]	687
[9098]	688	CALL lbc_lnk_multi( fr_iu , 'U', 1., fr_iv , 'V', 1. )
[2874]	689
[9019]	690	! set the snow+ice mass
	691	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	692	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	693	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	694	snwice_mass_b(:,:) = snwice_mass(:,:)
	695	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
[3193]	696	!
[2874]	697	END SUBROUTINE cice_sbc_out
	698
[3152]	699
[2874]	700	SUBROUTINE cice_sbc_hadgam( kt )
	701	!!---------------------------------------------------------------------
	702	!! * ROUTINE cice_sbc_hadgam *
[3040]	703	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
[2874]	704	!!
	705	!!
[9124]	706	!!---------------------------------------------------------------------
[2874]	707	INTEGER, INTENT( in ) :: kt ! ocean time step
[9124]	708	!!
[3625]	709	INTEGER :: jl ! dummy loop index
[3193]	710	INTEGER :: ierror
[9124]	711	!!---------------------------------------------------------------------
[9500]	712	#if defined key_asminc
	713	! Import fresh water and salt flux due to seaice da
	714	CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
	715	CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
	716	#endif
	717
[9498]	718	!
[2874]	719	! ! =========================== !
	720	! ! Prepare Coupling fields !
	721	! ! =========================== !
[9019]	722	!
	723	! x and y comp of ice velocity
	724	!
[3193]	725	CALL cice2nemo(uvel,u_ice,'F', -1. )
	726	CALL cice2nemo(vvel,v_ice,'F', -1. )
	727	!
[9019]	728	! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
	729	!
	730	! Snow and ice thicknesses (CO_2 and CO_3)
	731	!
	732	DO jl = 1, ncat
	733	CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
	734	CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
	735	END DO
[9500]	736
	737	#if ! defined key_cice4
	738	! Meltpond fraction and depth
	739	DO jl = 1,ncat
	740	CALL cice2nemo(apeffn(:,:,jl,:),a_p(:,:,jl),'T', 1. )
	741	CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),ht_p(:,:,jl),'T', 1. )
	742	ENDDO
	743	#endif
	744
	745
	746	! If using multilayers thermodynamics in CICE then get top layer temperature
	747	! and effective conductivity
	748	!! When using NEMO with CICE, this change requires use of
	749	!! one of the following two CICE branches:
	750	!! - at CICE5.0, hadax/r1015_GSI8_with_GSI7
	751	!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
	752	IF (heat_capacity) THEN
	753	DO jl = 1,ncat
	754	CALL cice2nemo(Tn_top(:,:,jl,:),tn_ice(:,:,jl),'T', 1. )
	755	CALL cice2nemo(keffn_top(:,:,jl,:),kn_ice(:,:,jl),'T', 1. )
	756	ENDDO
	757	! Convert surface temperature to Kelvin
	758	tn_ice(:,:,:)=tn_ice(:,:,:)+rt0
	759	ELSE
	760	tn_ice(:,:,:) = 0.0
	761	kn_ice(:,:,:) = 0.0
	762	ENDIF
	763
[9019]	764	!
[2874]	765	END SUBROUTINE cice_sbc_hadgam
	766
	767
	768	SUBROUTINE cice_sbc_final
	769	!!---------------------------------------------------------------------
	770	!! * ROUTINE cice_sbc_final *
	771	!! ** Purpose: Finalize CICE
	772	!!---------------------------------------------------------------------
[9124]	773	!
[2874]	774	IF(lwp) WRITE(numout,*)'cice_sbc_final'
[9124]	775	!
[3193]	776	CALL CICE_Finalize
[9124]	777	!
[2874]	778	END SUBROUTINE cice_sbc_final
	779
[9124]	780
[2874]	781	SUBROUTINE cice_sbc_force (kt)
	782	!!---------------------------------------------------------------------
	783	!! * ROUTINE cice_sbc_force *
	784	!! ** Purpose : Provide CICE forcing from files
	785	!!
	786	!!---------------------------------------------------------------------
	787	!! ** Method : READ monthly flux file in NetCDF files
	788	!!
	789	!! snowfall
	790	!! rainfall
	791	!! sublimation rate
	792	!! topmelt (category)
	793	!! botmelt (category)
	794	!!
	795	!! History :
	796	!!----------------------------------------------------------------------
	797	USE iom
[9124]	798	!!
	799	INTEGER, INTENT( in ) :: kt ! ocean time step
	800	!!
[2874]	801	INTEGER :: ierror ! return error code
	802	INTEGER :: ifpr ! dummy loop index
	803	!!
	804	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
	805	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
	806	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
	807	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
	808	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
	809	!!
	810	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
	811	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
	812	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
[4230]	813	INTEGER :: ios
[2874]	814	!!---------------------------------------------------------------------
	815
	816	! ! ====================== !
	817	IF( kt == nit000 ) THEN ! First call kt=nit000 !
	818	! ! ====================== !
[4990]	819	! namsbc_cice is not yet in the reference namelist
	820	! set file information (default values)
	821	cn_dir = './' ! directory in which the model is executed
	822
	823	! (NB: frequency positive => hours, negative => months)
	824	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask
	825	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file
	826	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	827	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' )
	828	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' )
	829	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	830	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	831	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	832	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	833	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	834	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
	835	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
	836	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
	837	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
	838	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
	839
[4230]	840	REWIND( numnam_ref ) ! Namelist namsbc_cice in reference namelist :
	841	READ ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
[9168]	842	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
[2874]	843
[4230]	844	REWIND( numnam_cfg ) ! Namelist namsbc_cice in configuration namelist : Parameters of the run
	845	READ ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
[9168]	846	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
[4624]	847	IF(lwm) WRITE ( numond, namsbc_cice )
[2874]	848
	849	! store namelist information in an array
	850	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
	851	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
	852	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
	853	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
	854	slf_i(jp_bot5) = sn_bot5
	855
	856	! set sf structure
	857	ALLOCATE( sf(jpfld), STAT=ierror )
	858	IF( ierror > 0 ) THEN
	859	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
	860	ENDIF
	861
	862	DO ifpr= 1, jpfld
	863	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
	864	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
	865	END DO
	866
	867	! fill sf with slf_i and control print
	868	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
	869	!
	870	ENDIF
	871
	872	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
	873	! ! input fields at the current time-step
	874
	875	! set the fluxes from read fields
	876	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
	877	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
[3040]	878	! May be better to do this conversion somewhere else
[2874]	879	qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
	880	topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
	881	topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
	882	topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
	883	topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
	884	topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
	885	botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
	886	botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
	887	botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
	888	botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
	889	botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
	890
	891	! control print (if less than 100 time-step asked)
	892	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
	893	WRITE(numout,*)
	894	WRITE(numout,*) ' read forcing fluxes for CICE OK'
	895	CALL FLUSH(numout)
	896	ENDIF
	897
	898	END SUBROUTINE cice_sbc_force
	899
	900	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
	901	!!---------------------------------------------------------------------
	902	!! * ROUTINE nemo2cice *
	903	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
	904	#if defined key_nemocice_decomp
	905	!!
	906	!! NEMO and CICE PE sub domains are identical, hence
	907	!! there is no need to gather or scatter data from
	908	!! one PE configuration to another.
	909	#else
	910	!! Automatically gather/scatter between
	911	!! different processors and blocks
	912	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
	913	!! B. Gather pn into global array (png)
	914	!! C. Map png into CICE global array (pcg)
	915	!! D. Scatter pcg to CICE blocks (pc) + update haloes
	916	#endif
	917	!!---------------------------------------------------------------------
[3193]	918	CHARACTER(len=1), INTENT( in ) :: &
	919	cd_type ! nature of pn grid-point
	920	! ! = T or F gridpoints
	921	REAL(wp), INTENT( in ) :: &
	922	psgn ! control of the sign change
	923	! ! =-1 , the sign is modified following the type of b.c. used
	924	! ! = 1 , no sign change
	925	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	926	#if !defined key_nemocice_decomp
[3625]	927	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
[3193]	928	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	929	#endif
[3193]	930	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
	931	INTEGER (int_kind) :: &
	932	field_type, &! id for type of field (scalar, vector, angle)
	933	grid_loc ! id for location on horizontal grid
[2874]	934	! (center, NEcorner, Nface, Eface)
	935
[3193]	936	INTEGER :: ji, jj, jn ! dummy loop indices
[5836]	937	!!---------------------------------------------------------------------
[2874]	938
[3193]	939	! A. Ensure all haloes are filled in NEMO field (pn)
[2874]	940
[3193]	941	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	942
	943	#if defined key_nemocice_decomp
	944
[3193]	945	! Copy local domain data from NEMO to CICE field
	946	pc(:,:,1)=0.0
[3625]	947	DO jj=2,ny_block-1
	948	DO ji=2,nx_block-1
	949	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
[3193]	950	ENDDO
	951	ENDDO
[2874]	952
	953	#else
	954
[3193]	955	! B. Gather pn into global array (png)
[2874]	956
[3193]	957	IF ( jpnij > 1) THEN
	958	CALL mppsync
	959	CALL mppgather (pn,0,png)
	960	CALL mppsync
	961	ELSE
	962	png(:,:,1)=pn(:,:)
	963	ENDIF
[2874]	964
[3193]	965	! C. Map png into CICE global array (pcg)
[2874]	966
	967	! Need to make sure this is robust to changes in NEMO halo rows....
	968	! (may be OK but not 100% sure)
	969
[3193]	970	IF (nproc==0) THEN
[2874]	971	! pcg(:,:)=0.0
[3193]	972	DO jn=1,jpnij
[3625]	973	DO jj=nldjt(jn),nlejt(jn)
	974	DO ji=nldit(jn),nleit(jn)
	975	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
[3193]	976	ENDDO
	977	ENDDO
	978	ENDDO
[3625]	979	DO jj=1,ny_global
	980	DO ji=1,nx_global
	981	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
	982	ENDDO
	983	ENDDO
[3193]	984	ENDIF
[2874]	985
	986	#endif
	987
[3193]	988	SELECT CASE ( cd_type )
	989	CASE ( 'T' )
	990	grid_loc=field_loc_center
	991	CASE ( 'F' )
	992	grid_loc=field_loc_NEcorner
	993	END SELECT
[2874]	994
[3193]	995	SELECT CASE ( NINT(psgn) )
	996	CASE ( -1 )
	997	field_type=field_type_vector
	998	CASE ( 1 )
	999	field_type=field_type_scalar
	1000	END SELECT
[2874]	1001
	1002	#if defined key_nemocice_decomp
[3193]	1003	! Ensure CICE halos are up to date
	1004	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1005	#else
[3193]	1006	! D. Scatter pcg to CICE blocks (pc) + update halos
	1007	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
[2874]	1008	#endif
	1009
	1010	END SUBROUTINE nemo2cice
	1011
	1012	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
	1013	!!---------------------------------------------------------------------
	1014	!! * ROUTINE cice2nemo *
	1015	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
	1016	#if defined key_nemocice_decomp
	1017	!!
	1018	!! NEMO and CICE PE sub domains are identical, hence
	1019	!! there is no need to gather or scatter data from
	1020	!! one PE configuration to another.
	1021	#else
	1022	!! Automatically deal with scatter/gather between
	1023	!! different processors and blocks
	1024	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
	1025	!! B. Map pcg into NEMO global array (png)
	1026	!! C. Scatter png into NEMO field (pn) for each processor
	1027	!! D. Ensure all haloes are filled in pn
	1028	#endif
	1029	!!---------------------------------------------------------------------
	1030
[3193]	1031	CHARACTER(len=1), INTENT( in ) :: &
	1032	cd_type ! nature of pn grid-point
	1033	! ! = T or F gridpoints
	1034	REAL(wp), INTENT( in ) :: &
	1035	psgn ! control of the sign change
	1036	! ! =-1 , the sign is modified following the type of b.c. used
	1037	! ! = 1 , no sign change
	1038	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	1039
	1040	#if defined key_nemocice_decomp
[3193]	1041	INTEGER (int_kind) :: &
	1042	field_type, & ! id for type of field (scalar, vector, angle)
	1043	grid_loc ! id for location on horizontal grid
	1044	! (center, NEcorner, Nface, Eface)
[2874]	1045	#else
[3193]	1046	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	1047	#endif
	1048
[3193]	1049	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
[2874]	1050
[3193]	1051	INTEGER :: ji, jj, jn ! dummy loop indices
[2874]	1052
	1053
	1054	#if defined key_nemocice_decomp
	1055
[3193]	1056	SELECT CASE ( cd_type )
	1057	CASE ( 'T' )
	1058	grid_loc=field_loc_center
	1059	CASE ( 'F' )
	1060	grid_loc=field_loc_NEcorner
	1061	END SELECT
[2874]	1062
[3193]	1063	SELECT CASE ( NINT(psgn) )
	1064	CASE ( -1 )
	1065	field_type=field_type_vector
	1066	CASE ( 1 )
	1067	field_type=field_type_scalar
	1068	END SELECT
[2874]	1069
[3193]	1070	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	1071
	1072
[3193]	1073	pn(:,:)=0.0
	1074	DO jj=1,jpjm1
	1075	DO ji=1,jpim1
[3625]	1076	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
[3193]	1077	ENDDO
	1078	ENDDO
[2874]	1079
	1080	#else
	1081
[3193]	1082	! A. Gather CICE blocks (pc) into global array (pcg)
[2874]	1083
[3193]	1084	CALL gather_global(pcg, pc, 0, distrb_info)
[2874]	1085
	1086	! B. Map pcg into NEMO global array (png)
	1087
	1088	! Need to make sure this is robust to changes in NEMO halo rows....
	1089	! (may be OK but not spent much time thinking about it)
[3625]	1090	! Note that non-existent pcg elements may be used below, but
	1091	! the lbclnk call on pn will replace these with sensible values
[2874]	1092
[3193]	1093	IF (nproc==0) THEN
	1094	png(:,:,:)=0.0
	1095	DO jn=1,jpnij
[3625]	1096	DO jj=nldjt(jn),nlejt(jn)
	1097	DO ji=nldit(jn),nleit(jn)
	1098	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
[3193]	1099	ENDDO
	1100	ENDDO
	1101	ENDDO
	1102	ENDIF
[2874]	1103
[3193]	1104	! C. Scatter png into NEMO field (pn) for each processor
[2874]	1105
[3193]	1106	IF ( jpnij > 1) THEN
	1107	CALL mppsync
	1108	CALL mppscatter (png,0,pn)
	1109	CALL mppsync
	1110	ELSE
	1111	pn(:,:)=png(:,:,1)
	1112	ENDIF
[2874]	1113
	1114	#endif
	1115
[3193]	1116	! D. Ensure all haloes are filled in pn
[2874]	1117
[3193]	1118	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	1119
	1120	END SUBROUTINE cice2nemo
	1121
	1122	#else
	1123	!!----------------------------------------------------------------------
	1124	!! Default option Dummy module NO CICE sea-ice model
	1125	!!----------------------------------------------------------------------
	1126	CONTAINS
	1127
[4990]	1128	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
[2874]	1129	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
	1130	END SUBROUTINE sbc_ice_cice
	1131
[4990]	1132	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
[2874]	1133	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
	1134	END SUBROUTINE cice_sbc_init
	1135
	1136	SUBROUTINE cice_sbc_final ! Dummy routine
	1137	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
	1138	END SUBROUTINE cice_sbc_final
	1139
	1140	#endif
	1141
	1142	!!======================================================================
	1143	END MODULE sbcice_cice

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

Download in other formats: