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

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

Merge in r8183 version of this branch (dev_r8183_GC_couple_pkg [8730:8734])

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