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_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 6675

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

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

Download in other formats: