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

source: branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 6041

Last change on this file since 6041 was 6041, checked in by timgraham, 8 years ago

Merged head of trunk into branch

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