source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9098

Last change on this file since 9098 was 9098, checked in by flavoni, 4 years ago

change lbc_lnk in lbc_lnk_multi

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