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 @ 12546

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

Adding precision specification in hardcoded reals and other modifications to allow compilation without forcing reals without precision specification to a certain value through compiler flags

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