source: branches/UKMO/dev_r5377_UKMO1_CICE_coupling_GSI7/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5522

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

Move the location where sstfrz is populated to after it is allocated bounds in sbc_ice_alloc.

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