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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8190

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

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

Download in other formats: