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

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

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

Context Navigation

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

Last change on this file since 5646 was 5646, checked in by dancopsey, 9 years ago
Moved sbc_cpl_alloc to be before sbc_cpl_init so that xcplmask in allocated dimensions before it is filled with data.
File size: 44.5 KB

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

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

Download in other formats: