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 branches/UKMO/r5936_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/r5936_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 7134

Last change on this file since 7134 was 7134, checked in by jcastill, 8 years ago

Version as in UKMO/dev_r5107_hadgem3_cplseq@5646

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