New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcice_cice.F90 in branches/UKMO/dev_r6912_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r6912_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9132

Last change on this file since 9132 was 9132, checked in by andmirek, 4 years ago

#1868 changes enabling coupling

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