Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sbcice_cice.F90 @ 11939

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_sit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11939

Download in other formats: