source: branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8478

Last change on this file since 8478 was 8478, checked in by jamrae, 3 years ago

Added call to cice2nemo for roughness length in cice_sbc_init.

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