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

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

Last change on this file was 9828, checked in by dancopsey, 6 years ago

Stop multiplying zevap_ice by -1.

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