source: branches/UKMO/dev_merge_2017_GC_couple_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9723

Last change on this file since 9723 was 9723, checked in by dancopsey, 2 years ago

Removed call top cice_sbc_hadgam is this is being done in the dev_merge_2017_CICE_interface branch. This means that we must always use the dev_merge_2017_CICE_interface branch with this branch.

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