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/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5479

Last change on this file since 5479 was 5479, checked in by cguiavarch, 9 years ago

Changes for sequence of coupling calls in HadGEM3.

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