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 NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcice_cice.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh

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