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_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/SBC – NEMO

source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcice_cice.F90 @ 13257

Last change on this file since 13257 was 13257, checked in by orioltp, 4 years ago

Updated with trunk at r13245 and small change allocating variables in icb_oce.F90.

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