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

Context Navigation

source: branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5646

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

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

Download in other formats: