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

Context Navigation

source: branches/UKMO/dev_r5518_med_test/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 6200

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

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

Download in other formats: