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

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 4147

Last change on this file since 4147 was 4147, checked in by cetlod, 10 years ago

merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169

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