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

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

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