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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 10276

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

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

Download in other formats: