New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcice_cice.F90 in branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 3261

Last change on this file since 3261 was 3261, checked in by charris, 12 years ago

Fix for restartability when running with CICE and key_vvl.

File size: 37.4 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 phycst, only : rcp, rau0
17   USE in_out_manager  ! I/O manager
18   USE lib_mpp         ! distributed memory computing library
19   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
20   USE wrk_nemo        ! work arrays
21   USE timing          ! Timing
22   USE daymod          ! calendar
23   USE fldread         ! read input fields
24
25   USE sbc_oce         ! Surface boundary condition: ocean fields
26   USE sbc_ice         ! Surface boundary condition: ice   fields
27   USE sbcblk_core     ! Surface boundary condition: CORE bulk
28   USE sbccpl
29
30   USE ice_kinds_mod
31   USE ice_blocks
32   USE ice_domain
33   USE ice_domain_size
34   USE ice_boundary
35   USE ice_constants
36   USE ice_gather_scatter
37   USE ice_calendar, only: dt
38   USE ice_state, only: aice,aicen,uvel,vvel,vsnon,vicen
39   USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
40                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm,     &
41                fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt,          &
42                flatn_f,fsurfn_f,fcondtopn_f,                    &
43                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
44                swvdr,swvdf,swidr,swidf
45   USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
46   USE ice_atmo, only: calc_strair
47   USE ice_therm_vertical, only: calc_Tsfc
48
49   USE CICE_InitMod
50   USE CICE_RunMod
51   USE CICE_FinalMod
52
53   IMPLICIT NONE
54   PRIVATE
55
56   !! * Routine accessibility
57   PUBLIC cice_sbc_init   ! routine called by sbc_init
58   PUBLIC cice_sbc_final  ! routine called by sbc_final
59   PUBLIC sbc_ice_cice    ! routine called by sbc
60
61   INTEGER , PARAMETER ::   jpfld   = 13   ! maximum number of files to read
62   INTEGER , PARAMETER ::   jp_snow = 1    ! index of snow file
63   INTEGER , PARAMETER ::   jp_rain = 2    ! index of rain file
64   INTEGER , PARAMETER ::   jp_sblm = 3    ! index of sublimation file
65   INTEGER , PARAMETER ::   jp_top1 = 4    ! index of category 1 topmelt file
66   INTEGER , PARAMETER ::   jp_top2 = 5    ! index of category 2 topmelt file
67   INTEGER , PARAMETER ::   jp_top3 = 6    ! index of category 3 topmelt file
68   INTEGER , PARAMETER ::   jp_top4 = 7    ! index of category 4 topmelt file
69   INTEGER , PARAMETER ::   jp_top5 = 8    ! index of category 5 topmelt file
70   INTEGER , PARAMETER ::   jp_bot1 = 9    ! index of category 1 botmelt file
71   INTEGER , PARAMETER ::   jp_bot2 = 10   ! index of category 2 botmelt file
72   INTEGER , PARAMETER ::   jp_bot3 = 11   ! index of category 3 botmelt file
73   INTEGER , PARAMETER ::   jp_bot4 = 12   ! index of category 4 botmelt file
74   INTEGER , PARAMETER ::   jp_bot5 = 13   ! index of category 5 botmelt file
75   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
76
77   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE ::   png     ! local array used in sbc_cice_ice
78
79   !! * Substitutions
80#  include "domzgr_substitute.h90"
81
82CONTAINS
83
84   INTEGER FUNCTION sbc_ice_cice_alloc()
85      !!----------------------------------------------------------------------
86      !!                ***  FUNCTION sbc_ice_cice_alloc  ***
87      !!----------------------------------------------------------------------
88      ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
89      IF( lk_mpp                 )   CALL mpp_sum ( sbc_ice_cice_alloc )
90      IF( sbc_ice_cice_alloc > 0 )   CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
91   END FUNCTION sbc_ice_cice_alloc
92
93   SUBROUTINE sbc_ice_cice( kt, nsbc )
94      !!---------------------------------------------------------------------
95      !!                  ***  ROUTINE sbc_ice_cice  ***
96      !!                   
97      !! ** Purpose :   update the ocean surface boundary condition via the
98      !!                CICE Sea Ice Model time stepping
99      !!
100      !! ** Method  : - Get any extra forcing fields for CICE 
101      !!              - Prepare forcing fields
102      !!              - CICE model time stepping
103      !!              - call the routine that computes mass and
104      !!                heat fluxes at the ice/ocean interface
105      !!
106      !! ** Action  : - time evolution of the CICE sea-ice model
107      !!              - update all sbc variables below sea-ice:
108      !!                utau, vtau, qns , qsr, emp , emps
109      !!---------------------------------------------------------------------
110      INTEGER, INTENT(in) ::   kt      ! ocean time step
111      INTEGER, INTENT(in) ::   nsbc    ! surface forcing type
112      !!----------------------------------------------------------------------
113      !
114      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_cice')
115      !
116      !                                        !----------------------!
117      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
118         !                                     !----------------------!
119
120         ! Make sure any fluxes required for CICE are set
121         IF ( nsbc == 2 )  THEN
122            CALL cice_sbc_force(kt)
123         ELSE IF ( nsbc == 5 ) THEN
124            CALL sbc_cpl_ice_flx( 1.0-fr_i  )
125         ENDIF
126
127         CALL cice_sbc_in ( kt, nsbc )
128         CALL CICE_Run
129         CALL cice_sbc_out ( kt, nsbc )
130
131         IF ( nsbc == 5 )  CALL cice_sbc_hadgam(kt+1)
132
133      ENDIF                                          ! End sea-ice time step only
134      !
135      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_cice')
136
137   END SUBROUTINE sbc_ice_cice
138
139   SUBROUTINE cice_sbc_init (nsbc)
140      !!---------------------------------------------------------------------
141      !!                    ***  ROUTINE cice_sbc_init  ***
142      !! ** Purpose: Initialise ice related fields for NEMO and coupling
143      !!
144      INTEGER, INTENT( in  ) ::   nsbc ! surface forcing type
145      !!---------------------------------------------------------------------
146
147      INTEGER  ::   ji, jj, jpl                        ! dummy loop indices
148
149      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_init')
150      !
151      IF(lwp) WRITE(numout,*)'cice_sbc_init'
152
153! Initialize CICE
154      CALL CICE_Initialize
155
156! Do some CICE consistency checks
157      IF ( (nsbc == 2) .OR. (nsbc == 5) ) THEN
158         IF ( calc_strair .OR. calc_Tsfc ) THEN
159            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
160         ENDIF
161      ELSEIF (nsbc == 4) THEN
162         IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
163            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
164         ENDIF
165      ENDIF
166
167
168! allocate sbc_ice and sbc_cice arrays
169      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' )
170      IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
171
172! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
173      IF( .NOT. ln_rstart ) THEN
174         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
175         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
176      ENDIF
177
178      fr_iu(:,:)=0.0
179      fr_iv(:,:)=0.0
180
181      CALL cice2nemo(aice,fr_i, 'T', 1. )
182      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
183         DO jpl=1,ncat
184            CALL cice2nemo(aicen(:,:,jpl,:),a_i(:,:,jpl), 'T', 1. )
185         ENDDO
186      ENDIF
187
188! T point to U point
189! T point to V point
190      DO jj=1,jpjm1
191         DO ji=1,jpim1
192            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
193            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
194         ENDDO
195      ENDDO
196
197      CALL lbc_lnk ( fr_iu , 'U', 1. )
198      CALL lbc_lnk ( fr_iv , 'V', 1. )
199      !
200      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init')
201      !
202   END SUBROUTINE cice_sbc_init
203
204   
205   SUBROUTINE cice_sbc_in (kt, nsbc)
206      !!---------------------------------------------------------------------
207      !!                    ***  ROUTINE cice_sbc_in  ***
208      !! ** Purpose: Set coupling fields and pass to CICE
209      !!---------------------------------------------------------------------
210      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
211      INTEGER, INTENT(in   ) ::   nsbc ! surface forcing type
212
213      INTEGER  ::   ji, jj, jpl                   ! dummy loop indices     
214      REAL(wp), DIMENSION(:,:), POINTER :: ztmp
215      REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
216      !!---------------------------------------------------------------------
217
218      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_in')
219      !
220      CALL wrk_alloc( jpi,jpj, ztmp )
221      CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
222
223      IF( kt == nit000 )  THEN
224         IF(lwp) WRITE(numout,*)'cice_sbc_in'
225      ENDIF
226
227      ztmp(:,:)=0.0
228
229! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
230! the first time-step)
231
232! forced and coupled case
233
234      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
235
236         ztmpn(:,:,:)=0.0
237
238! x comp of wind stress (CI_1)
239! U point to F point
240         DO jj=1,jpjm1
241            DO ji=1,jpi
242               ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
243                                    + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
244            ENDDO
245         ENDDO
246         CALL nemo2cice(ztmp,strax,'F', -1. )
247
248! y comp of wind stress (CI_2)
249! V point to F point
250         DO jj=1,jpj
251            DO ji=1,jpim1
252               ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
253                                    + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
254            ENDDO
255         ENDDO
256         CALL nemo2cice(ztmp,stray,'F', -1. )
257
258! Surface downward latent heat flux (CI_5)
259         IF (nsbc == 2) THEN
260            DO jpl=1,ncat
261               ztmpn(:,:,jpl)=qla_ice(:,:,1)*a_i(:,:,jpl)
262            ENDDO
263         ELSE
264! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
265            qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub
266! End of temporary code
267            DO jj=1,jpj
268               DO ji=1,jpi
269                  IF (fr_i(ji,jj).eq.0.0) THEN
270                     DO jpl=1,ncat
271                        ztmpn(ji,jj,jpl)=0.0
272                     ENDDO
273                     ! This will then be conserved in CICE
274                     ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
275                  ELSE
276                     DO jpl=1,ncat
277                        ztmpn(ji,jj,jpl)=qla_ice(ji,jj,1)*a_i(ji,jj,jpl)/fr_i(ji,jj)
278                     ENDDO
279                  ENDIF
280               ENDDO
281            ENDDO
282         ENDIF
283         DO jpl=1,ncat
284            CALL nemo2cice(ztmpn(:,:,jpl),flatn_f(:,:,jpl,:),'T', 1. )
285
286! GBM conductive flux through ice (CI_6)
287!  Convert to GBM
288            IF (nsbc == 2) THEN
289               ztmp(:,:) = botmelt(:,:,jpl)*a_i(:,:,jpl)
290            ELSE
291               ztmp(:,:) = botmelt(:,:,jpl)
292            ENDIF
293            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jpl,:),'T', 1. )
294
295! GBM surface heat flux (CI_7)
296!  Convert to GBM
297            IF (nsbc == 2) THEN
298               ztmp(:,:) = (topmelt(:,:,jpl)+botmelt(:,:,jpl))*a_i(:,:,jpl) 
299            ELSE
300               ztmp(:,:) = (topmelt(:,:,jpl)+botmelt(:,:,jpl))
301            ENDIF
302            CALL nemo2cice(ztmp,fsurfn_f(:,:,jpl,:),'T', 1. )
303         ENDDO
304
305      ELSE IF (nsbc == 4) THEN
306
307! Pass CORE forcing fields to CICE (which will calculate heat fluxes etc itself)
308! x comp and y comp of atmosphere surface wind (CICE expects on T points)
309         ztmp(:,:) = wndi_ice(:,:)
310         CALL nemo2cice(ztmp,uatm,'T', -1. )
311         ztmp(:,:) = wndj_ice(:,:)
312         CALL nemo2cice(ztmp,vatm,'T', -1. )
313         ztmp(:,:) = SQRT ( wndi_ice(:,:)**2 + wndj_ice(:,:)**2 )
314         CALL nemo2cice(ztmp,wind,'T', 1. )    ! Wind speed (m/s)
315         ztmp(:,:) = qsr_ice(:,:,1)
316         CALL nemo2cice(ztmp,fsw,'T', 1. )     ! Incoming short-wave (W/m^2)
317         ztmp(:,:) = qlw_ice(:,:,1)
318         CALL nemo2cice(ztmp,flw,'T', 1. )     ! Incoming long-wave (W/m^2)
319         ztmp(:,:) = tatm_ice(:,:)
320         CALL nemo2cice(ztmp,Tair,'T', 1. )    ! Air temperature (K)
321         CALL nemo2cice(ztmp,potT,'T', 1. )    ! Potential temp (K)
322! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows 
323         ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )   
324                                               ! Constant (101000.) atm pressure assumed
325         CALL nemo2cice(ztmp,rhoa,'T', 1. )    ! Air density (kg/m^3)
326         ztmp(:,:) = qatm_ice(:,:)
327         CALL nemo2cice(ztmp,Qa,'T', 1. )      ! Specific humidity (kg/kg)
328         ztmp(:,:)=10.0
329         CALL nemo2cice(ztmp,zlvl,'T', 1. )    ! Atmos level height (m)
330
331! May want to check all values are physically realistic (as in CICE routine
332! prepare_forcing)?
333
334! Divide shortwave into spectral bands (as in prepare_forcing)
335         ztmp(:,:)=qsr_ice(:,:,1)*frcvdr       ! visible direct
336         CALL nemo2cice(ztmp,swvdr,'T', 1. )             
337         ztmp(:,:)=qsr_ice(:,:,1)*frcvdf       ! visible diffuse
338         CALL nemo2cice(ztmp,swvdf,'T', 1. )             
339         ztmp(:,:)=qsr_ice(:,:,1)*frcidr       ! near IR direct
340         CALL nemo2cice(ztmp,swidr,'T', 1. )
341         ztmp(:,:)=qsr_ice(:,:,1)*frcidf       ! near IR diffuse
342         CALL nemo2cice(ztmp,swidf,'T', 1. )
343
344      ENDIF
345
346! Snowfall
347! Ensure fsnow is positive (as in CICE routine prepare_forcing) 
348      ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0) 
349      CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
350
351! Rainfall
352      ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
353      CALL nemo2cice(ztmp,frain,'T', 1. ) 
354
355! Freezing/melting potential
356! Calculated over NEMO leapfrog timestep (hence 2*dt), so uses tb in
357! frzmlt which is then applied in going from tb to ta.
358! May be better using sst_m if not coupling to CICE every time-step
359
360!      nfrzmlt(:,:)=rau0*rcp*fse3t(:,:,1)*(Tocnfrz-sst_m(:,:))/(2.0*dt)
361      nfrzmlt(:,:)=rau0*rcp*fse3t_b(:,:,1)*(Tocnfrz-tsb(:,:,1,jp_tem))/(2.0*dt)
362
363      ztmp(:,:) = nfrzmlt(:,:)
364      CALL nemo2cice(ztmp,frzmlt,'T', 1. )
365
366! SST  and SSS
367
368      CALL nemo2cice(sst_m,sst,'T', 1. )
369      CALL nemo2cice(sss_m,sss,'T', 1. )
370
371! x comp and y comp of surface ocean current
372! U point to F point
373      DO jj=1,jpjm1
374         DO ji=1,jpi
375            ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
376         ENDDO
377      ENDDO
378      CALL nemo2cice(ztmp,uocn,'F', -1. )
379
380! V point to F point
381      DO jj=1,jpj
382         DO ji=1,jpim1
383            ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
384         ENDDO
385      ENDDO
386      CALL nemo2cice(ztmp,vocn,'F', -1. )
387
388! x comp and y comp of sea surface slope (on F points)
389! T point to F point
390      DO jj=1,jpjm1
391         DO ji=1,jpim1
392            ztmp(ji,jj)=0.5 * (  (ssh_m(ji+1,jj  )-ssh_m(ji,jj  ))/e1u(ji,jj  )   &
393                               + (ssh_m(ji+1,jj+1)-ssh_m(ji,jj+1))/e1u(ji,jj+1) ) & 
394                            *  fmask(ji,jj,1)
395         ENDDO
396      ENDDO
397      CALL nemo2cice(ztmp,ss_tltx,'F', -1. )
398
399! T point to F point
400      DO jj=1,jpjm1
401         DO ji=1,jpim1
402            ztmp(ji,jj)=0.5 * (  (ssh_m(ji  ,jj+1)-ssh_m(ji  ,jj))/e2v(ji  ,jj)   &
403                               + (ssh_m(ji+1,jj+1)-ssh_m(ji+1,jj))/e2v(ji+1,jj) ) &
404                            *  fmask(ji,jj,1)
405         ENDDO
406      ENDDO
407      CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
408
409      CALL wrk_dealloc( jpi,jpj, ztmp )
410      CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
411      !
412      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_in')
413      !
414   END SUBROUTINE cice_sbc_in
415
416
417   SUBROUTINE cice_sbc_out (kt,nsbc)
418      !!---------------------------------------------------------------------
419      !!                    ***  ROUTINE cice_sbc_out  ***
420      !! ** Purpose: Get fields from CICE and set surface fields for NEMO
421      !!---------------------------------------------------------------------
422      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
423      INTEGER, INTENT( in  ) ::   nsbc ! surface forcing type
424     
425      INTEGER  ::   ji, jj, jpl                 ! dummy loop indices
426      REAL(wp), DIMENSION(:,:), POINTER :: ztmp
427      !!---------------------------------------------------------------------
428
429      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_out')
430      !
431      CALL wrk_alloc( jpi,jpj, ztmp )
432     
433      IF( kt == nit000 )  THEN
434         IF(lwp) WRITE(numout,*)'cice_sbc_out'
435      ENDIF
436     
437! x comp of ocean-ice stress
438      CALL cice2nemo(strocnx,ztmp,'F', -1. )
439      ss_iou(:,:)=0.0
440! F point to U point
441      DO jj=2,jpjm1
442         DO ji=2,jpim1
443            ss_iou(ji,jj) = 0.5 * ( ztmp(ji,jj-1) + ztmp(ji,jj) ) * umask(ji,jj,1)
444         ENDDO
445      ENDDO
446      CALL lbc_lnk( ss_iou , 'U', -1. )
447
448! y comp of ocean-ice stress
449      CALL cice2nemo(strocny,ztmp,'F', -1. )
450      ss_iov(:,:)=0.0
451! F point to V point
452
453      DO jj=1,jpjm1
454         DO ji=2,jpim1
455            ss_iov(ji,jj) = 0.5 * ( ztmp(ji-1,jj) + ztmp(ji,jj) ) * vmask(ji,jj,1)
456         ENDDO
457      ENDDO
458      CALL lbc_lnk( ss_iov , 'V', -1. )
459
460! x and y comps of surface stress
461! Combine wind stress and ocean-ice stress
462! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
463
464      utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
465      vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)     
466
467! Freshwater fluxes
468
469      IF (nsbc == 2) THEN
470! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
471! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below
472! Not ideal since aice won't be the same as in the atmosphere. 
473! Better to use evap and tprecip? (but for now don't read in evap in this case)
474         emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
475      ELSE IF (nsbc == 4) THEN
476         emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)       
477      ELSE IF (nsbc ==5) THEN
478! emp_tot is set in sbc_cpl_ice_flx (call from cice_sbc_in above)
479         emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:) 
480      ENDIF
481
482! Subtract fluxes from CICE to get freshwater equivalent flux used in
483! salinity calculation
484      CALL cice2nemo(fresh_gbm,ztmp,'T', 1. )
485      emps(:,:)=emp(:,:)-ztmp(:,:)
486! Note the 1000.0 is to convert from kg salt to g salt (needed for PSU)
487      CALL cice2nemo(fsalt_gbm,ztmp,'T', 1. )
488      DO jj=1,jpj
489         DO ji=1,jpi
490            IF (sss_m(ji,jj).gt.0.0) THEN
491               emps(ji,jj)=emps(ji,jj)+ztmp(ji,jj)*1000.0/sss_m(ji,jj)
492            ENDIF
493         ENDDO
494      ENDDO
495
496! No longer remove precip over ice from free surface calculation on basis that the
497! weight of the precip will affect the free surface even if it falls on the ice
498! (same to the argument that freezing / melting of ice doesn't change the free surface)
499! Sublimation from the ice is treated in a similar way (included in emp but not emps) 
500!
501! This will need re-thinking in the variable volume case and if ice is 'embedded' in the
502! ocean rather than floating on top
503
504      emp(:,:)  = emp(:,:) - tprecip(:,:)*fr_i(:,:)
505
506! Take sublimation into account
507      IF (nsbc == 5 ) THEN
508         emp(:,:) = emp(:,:) + ( emp_ice(:,:) + sprecip(:,:) )
509      ELSE IF (nsbc == 2 ) THEN
510         emp(:,:) = emp(:,:) - qla_ice(:,:,1) / Lsub
511      ENDIF
512
513      CALL lbc_lnk( emp , 'T', 1. )
514      CALL lbc_lnk( emps , '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.