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

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5516

Last change on this file since 5516 was 5516, checked in by smasson, 9 years ago

work_alloc/work_dealloc errors

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