New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcice_cice.F90 in NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcice_cice.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 3 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496

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