source: branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 6458

Last change on this file since 6458 was 6458, checked in by dancopsey, 5 years ago

Mame eos_fzp be called as a subroutine instead of a function to be compatible with changes in nemo_v3_6_STABLE_copy revision 6436

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