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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcice_cice.F90 in branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8384

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

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

Download in other formats: