source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 3402

Last change on this file since 3402 was 3402, checked in by acc, 9 years ago

Branch: dev_r3385_NOCS04_HAMF; #665. Stage 2 of 2012 development: suppression of emps array and introduction of sfx (salt flux) array with associated code to setup the options for embedding the seaice into the ocean

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