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 @ 3152

Last change on this file since 3152 was 3152, checked in by smasson, 12 years ago

dev_NEMO_MERGE_2011: new dynamical allocation in IOM and SBC

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