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

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 3508

Last change on this file since 3508 was 3508, checked in by charris, 12 years ago

First attempt at changes for CICE embedded sea ice. The nn_ice_embd=0 option has been prevented (as for LIM3). Changes to sbcice_cice include some other fixes
and cosmetic changes as in dev_3352_UKMO8_CICE (care will be required merging with this branch).

File size: 41.9 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!         REWIND ( numnam )               ! ... at some point might read in from NEMO namelist?
762!         READ   ( numnam, namsbc_cice )
763
764         ! store namelist information in an array
765         slf_i(jp_snow) = sn_snow   ;   slf_i(jp_rain) = sn_rain   ;   slf_i(jp_sblm) = sn_sblm
766         slf_i(jp_top1) = sn_top1   ;   slf_i(jp_top2) = sn_top2   ;   slf_i(jp_top3) = sn_top3
767         slf_i(jp_top4) = sn_top4   ;   slf_i(jp_top5) = sn_top5   ;   slf_i(jp_bot1) = sn_bot1
768         slf_i(jp_bot2) = sn_bot2   ;   slf_i(jp_bot3) = sn_bot3   ;   slf_i(jp_bot4) = sn_bot4
769         slf_i(jp_bot5) = sn_bot5
770         
771         ! set sf structure
772         ALLOCATE( sf(jpfld), STAT=ierror )
773         IF( ierror > 0 ) THEN
774            CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' )   ;   RETURN
775         ENDIF
776
777         DO ifpr= 1, jpfld
778            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
779            ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
780         END DO
781
782         ! fill sf with slf_i and control print
783         CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
784         !
785      ENDIF
786
787      CALL fld_read( kt, nn_fsbc, sf )           ! Read input fields and provides the
788      !                                          ! input fields at the current time-step
789
790      ! set the fluxes from read fields
791      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
792      tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
793! May be better to do this conversion somewhere else
794      qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
795      topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
796      topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
797      topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
798      topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
799      topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
800      botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
801      botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
802      botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
803      botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
804      botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
805
806      ! control print (if less than 100 time-step asked)
807      IF( nitend-nit000 <= 100 .AND. lwp ) THEN
808         WRITE(numout,*) 
809         WRITE(numout,*) '        read forcing fluxes for CICE OK'
810         CALL FLUSH(numout)
811      ENDIF
812
813   END SUBROUTINE cice_sbc_force
814
815   SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
816      !!---------------------------------------------------------------------
817      !!                    ***  ROUTINE nemo2cice  ***
818      !! ** Purpose :   Transfer field in NEMO array to field in CICE array. 
819#if defined key_nemocice_decomp
820      !!             
821      !!                NEMO and CICE PE sub domains are identical, hence
822      !!                there is no need to gather or scatter data from
823      !!                one PE configuration to another.
824#else
825      !!                Automatically gather/scatter between
826      !!                different processors and blocks
827      !! ** Method :    A. Ensure all haloes are filled in NEMO field (pn)
828      !!                B. Gather pn into global array (png)
829      !!                C. Map png into CICE global array (pcg)
830      !!                D. Scatter pcg to CICE blocks (pc) + update haloes 
831#endif
832      !!---------------------------------------------------------------------
833
834      CHARACTER(len=1), INTENT( in ) ::   &
835          cd_type       ! nature of pn grid-point
836          !             !   = T or F gridpoints
837      REAL(wp), INTENT( in ) ::   &
838          psgn          ! control of the sign change
839          !             !   =-1 , the sign is modified following the type of b.c. used
840          !             !   = 1 , no sign change
841      REAL(wp), DIMENSION(jpi,jpj) :: pn
842#if !defined key_nemocice_decomp
843      REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
844      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
845#endif
846      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
847      INTEGER (int_kind) :: &
848         field_type,        &! id for type of field (scalar, vector, angle)
849         grid_loc            ! id for location on horizontal grid
850                            !  (center, NEcorner, Nface, Eface)
851
852      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
853
854!     A. Ensure all haloes are filled in NEMO field (pn)
855
856      CALL lbc_lnk( pn , cd_type, psgn )
857
858#if defined key_nemocice_decomp
859
860      ! Copy local domain data from NEMO to CICE field
861      pc(:,:,1)=0.0
862      DO jj=2,ny_block-1
863         DO ji=2,nx_block-1
864            pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
865         ENDDO
866      ENDDO
867
868#else
869
870!     B. Gather pn into global array (png)
871
872      IF ( jpnij > 1) THEN
873         CALL mppsync
874         CALL mppgather (pn,0,png) 
875         CALL mppsync
876      ELSE
877         png(:,:,1)=pn(:,:)
878      ENDIF
879
880!     C. Map png into CICE global array (pcg)
881
882! Need to make sure this is robust to changes in NEMO halo rows....
883! (may be OK but not 100% sure)
884
885      IF (nproc==0) THEN     
886!        pcg(:,:)=0.0
887         DO jn=1,jpnij
888            DO jj=nldjt(jn),nlejt(jn)
889               DO ji=nldit(jn),nleit(jn)
890                  png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
891               ENDDO
892            ENDDO
893         ENDDO
894         DO jj=1,ny_global
895            DO ji=1,nx_global
896               pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
897            ENDDO
898         ENDDO
899      ENDIF
900
901#endif
902
903      SELECT CASE ( cd_type )
904         CASE ( 'T' )
905            grid_loc=field_loc_center
906         CASE ( 'F' )                             
907            grid_loc=field_loc_NEcorner
908      END SELECT
909
910      SELECT CASE ( NINT(psgn) )
911         CASE ( -1 )
912            field_type=field_type_vector
913         CASE ( 1 )                             
914            field_type=field_type_scalar
915      END SELECT
916
917#if defined key_nemocice_decomp
918      ! Ensure CICE halos are up to date
919      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
920#else
921!     D. Scatter pcg to CICE blocks (pc) + update halos
922      CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
923#endif
924
925   END SUBROUTINE nemo2cice
926
927   SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
928      !!---------------------------------------------------------------------
929      !!                    ***  ROUTINE cice2nemo  ***
930      !! ** Purpose :   Transfer field in CICE array to field in NEMO array.
931#if defined key_nemocice_decomp
932      !!             
933      !!                NEMO and CICE PE sub domains are identical, hence
934      !!                there is no need to gather or scatter data from
935      !!                one PE configuration to another.
936#else 
937      !!                Automatically deal with scatter/gather between
938      !!                different processors and blocks
939      !! ** Method :    A. Gather CICE blocks (pc) into global array (pcg)
940      !!                B. Map pcg into NEMO global array (png)
941      !!                C. Scatter png into NEMO field (pn) for each processor
942      !!                D. Ensure all haloes are filled in pn
943#endif
944      !!---------------------------------------------------------------------
945
946      CHARACTER(len=1), INTENT( in ) ::   &
947          cd_type       ! nature of pn grid-point
948          !             !   = T or F gridpoints
949      REAL(wp), INTENT( in ) ::   &
950          psgn          ! control of the sign change
951          !             !   =-1 , the sign is modified following the type of b.c. used
952          !             !   = 1 , no sign change
953      REAL(wp), DIMENSION(jpi,jpj) :: pn
954
955#if defined key_nemocice_decomp
956      INTEGER (int_kind) :: &
957         field_type,        & ! id for type of field (scalar, vector, angle)
958         grid_loc             ! id for location on horizontal grid
959                              ! (center, NEcorner, Nface, Eface)
960#else
961      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
962#endif
963
964      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
965
966      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
967
968
969#if defined key_nemocice_decomp
970
971      SELECT CASE ( cd_type )
972         CASE ( 'T' )
973            grid_loc=field_loc_center
974         CASE ( 'F' )                             
975            grid_loc=field_loc_NEcorner
976      END SELECT
977
978      SELECT CASE ( NINT(psgn) )
979         CASE ( -1 )
980            field_type=field_type_vector
981         CASE ( 1 )                             
982            field_type=field_type_scalar
983      END SELECT
984
985      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
986
987
988      pn(:,:)=0.0
989      DO jj=1,jpjm1
990         DO ji=1,jpim1
991            pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
992         ENDDO
993      ENDDO
994
995#else
996
997!      A. Gather CICE blocks (pc) into global array (pcg)
998
999      CALL gather_global(pcg, pc, 0, distrb_info)
1000
1001!     B. Map pcg into NEMO global array (png)
1002
1003! Need to make sure this is robust to changes in NEMO halo rows....
1004! (may be OK but not spent much time thinking about it)
1005! Note that non-existent pcg elements may be used below, but
1006! the lbclnk call on pn will replace these with sensible values
1007
1008      IF (nproc==0) THEN
1009         png(:,:,:)=0.0
1010         DO jn=1,jpnij
1011            DO jj=nldjt(jn),nlejt(jn)
1012               DO ji=nldit(jn),nleit(jn)
1013                  png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
1014               ENDDO
1015            ENDDO
1016         ENDDO
1017      ENDIF
1018
1019!     C. Scatter png into NEMO field (pn) for each processor
1020
1021      IF ( jpnij > 1) THEN
1022         CALL mppsync
1023         CALL mppscatter (png,0,pn) 
1024         CALL mppsync
1025      ELSE
1026         pn(:,:)=png(:,:,1)
1027      ENDIF
1028
1029#endif
1030
1031!     D. Ensure all haloes are filled in pn
1032
1033      CALL lbc_lnk( pn , cd_type, psgn )
1034
1035   END SUBROUTINE cice2nemo
1036
1037#else
1038   !!----------------------------------------------------------------------
1039   !!   Default option           Dummy module         NO CICE sea-ice model
1040   !!----------------------------------------------------------------------
1041CONTAINS
1042
1043   SUBROUTINE sbc_ice_cice ( kt, nsbc )     ! Dummy routine
1044      WRITE(*,*) 'sbc_ice_cice: You should not have seen this print! error?', kt
1045   END SUBROUTINE sbc_ice_cice
1046
1047   SUBROUTINE cice_sbc_init (nsbc)    ! Dummy routine
1048      WRITE(*,*) 'cice_sbc_init: You should not have seen this print! error?'
1049   END SUBROUTINE cice_sbc_init
1050
1051   SUBROUTINE cice_sbc_final     ! Dummy routine
1052      WRITE(*,*) 'cice_sbc_final: You should not have seen this print! error?'
1053   END SUBROUTINE cice_sbc_final
1054
1055#endif
1056
1057   !!======================================================================
1058END MODULE sbcice_cice
Note: See TracBrowser for help on using the repository browser.