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

Context Navigation

source: branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5086

Last change on this file since 5086 was 5086, checked in by timgraham, 9 years ago

Merged head of trunk into branch in preparation for putting code back onto the trunk
In working copy ran the command:
svn merge svn+sshtimgraham@…/ipsl/forge/projets/nemo/svn/trunk

Also recompiled NEMO_book.pdf with merged input files

File size: 43.8 KB

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

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

Download in other formats: