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

source: branches/2011/dev_r2802_UKMO8_cice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 2884

Last change on this file since 2884 was 2884, checked in by charris, 13 years ago

#662 Tidying + changes to make sure NEMO-CICE model will compile OK.

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