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.
sbccpl.F90 in branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 4148

Last change on this file since 4148 was 4148, checked in by cetlod, 10 years ago

merge in trunk changes between r3853 and r3940 and commit the changes, see ticket #1169

  • Property svn:keywords set to Id
File size: 95.4 KB
Line 
1MODULE sbccpl
2   !!======================================================================
3   !!                       ***  MODULE  sbccpl  ***
4   !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode
5   !!======================================================================
6   !! History :  2.0  ! 2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
7   !!            3.0  ! 2008-02  (G. Madec, C Talandier)  surface module
8   !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
9   !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
10   !!----------------------------------------------------------------------
11#if defined key_oasis3 || defined key_oasis4
12   !!----------------------------------------------------------------------
13   !!   'key_oasis3' or 'key_oasis4'   Coupled Ocean/Atmosphere formulation
14   !!----------------------------------------------------------------------
15   !!   namsbc_cpl      : coupled formulation namlist
16   !!   sbc_cpl_init    : initialisation of the coupled exchanges
17   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
18   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
19   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
20   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
21   !!   sbc_cpl_snd     : send     fields to the atmosphere
22   !!----------------------------------------------------------------------
23   USE dom_oce         ! ocean space and time domain
24   USE sbc_oce         ! Surface boundary condition: ocean fields
25   USE sbc_ice         ! Surface boundary condition: ice fields
26   USE sbcdcy          ! surface boundary condition: diurnal cycle
27   USE phycst          ! physical constants
28#if defined key_lim3
29   USE par_ice         ! ice parameters
30   USE ice             ! ice variables
31#endif
32#if defined key_lim2
33   USE par_ice_2       ! ice parameters
34   USE ice_2           ! ice variables
35#endif
36#if defined key_oasis3
37   USE cpl_oasis3      ! OASIS3 coupling
38#endif
39#if defined key_oasis4
40   USE cpl_oasis4      ! OASIS4 coupling
41#endif
42   USE geo2ocean       !
43   USE oce   , ONLY : tsn, un, vn
44   USE albedo          !
45   USE in_out_manager  ! I/O manager
46   USE iom             ! NetCDF library
47   USE lib_mpp         ! distribued memory computing library
48   USE wrk_nemo        ! work arrays
49   USE timing          ! Timing
50   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
51#if defined key_cpl_carbon_cycle
52   USE p4zflx, ONLY : oce_co2
53#endif
54   USE diaar5, ONLY :   lk_diaar5
55#if defined key_cice
56   USE ice_domain_size, only: ncat
57#endif
58   IMPLICIT NONE
59   PRIVATE
60
61   PUBLIC   sbc_cpl_rcv        ! routine called by sbc_ice_lim(_2).F90
62   PUBLIC   sbc_cpl_snd        ! routine called by step.F90
63   PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
64   PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).F90
65
66   INTEGER, PARAMETER ::   jpr_otx1   =  1            ! 3 atmosphere-ocean stress components on grid 1
67   INTEGER, PARAMETER ::   jpr_oty1   =  2            !
68   INTEGER, PARAMETER ::   jpr_otz1   =  3            !
69   INTEGER, PARAMETER ::   jpr_otx2   =  4            ! 3 atmosphere-ocean stress components on grid 2
70   INTEGER, PARAMETER ::   jpr_oty2   =  5            !
71   INTEGER, PARAMETER ::   jpr_otz2   =  6            !
72   INTEGER, PARAMETER ::   jpr_itx1   =  7            ! 3 atmosphere-ice   stress components on grid 1
73   INTEGER, PARAMETER ::   jpr_ity1   =  8            !
74   INTEGER, PARAMETER ::   jpr_itz1   =  9            !
75   INTEGER, PARAMETER ::   jpr_itx2   = 10            ! 3 atmosphere-ice   stress components on grid 2
76   INTEGER, PARAMETER ::   jpr_ity2   = 11            !
77   INTEGER, PARAMETER ::   jpr_itz2   = 12            !
78   INTEGER, PARAMETER ::   jpr_qsroce = 13            ! Qsr above the ocean
79   INTEGER, PARAMETER ::   jpr_qsrice = 14            ! Qsr above the ice
80   INTEGER, PARAMETER ::   jpr_qsrmix = 15 
81   INTEGER, PARAMETER ::   jpr_qnsoce = 16            ! Qns above the ocean
82   INTEGER, PARAMETER ::   jpr_qnsice = 17            ! Qns above the ice
83   INTEGER, PARAMETER ::   jpr_qnsmix = 18
84   INTEGER, PARAMETER ::   jpr_rain   = 19            ! total liquid precipitation (rain)
85   INTEGER, PARAMETER ::   jpr_snow   = 20            ! solid precipitation over the ocean (snow)
86   INTEGER, PARAMETER ::   jpr_tevp   = 21            ! total evaporation
87   INTEGER, PARAMETER ::   jpr_ievp   = 22            ! solid evaporation (sublimation)
88   INTEGER, PARAMETER ::   jpr_sbpr   = 23            ! sublimation - liquid precipitation - solid precipitation
89   INTEGER, PARAMETER ::   jpr_semp   = 24            ! solid freshwater budget (sublimation - snow)
90   INTEGER, PARAMETER ::   jpr_oemp   = 25            ! ocean freshwater budget (evap - precip)
91   INTEGER, PARAMETER ::   jpr_w10m   = 26            ! 10m wind
92   INTEGER, PARAMETER ::   jpr_dqnsdt = 27            ! d(Q non solar)/d(temperature)
93   INTEGER, PARAMETER ::   jpr_rnf    = 28            ! runoffs
94   INTEGER, PARAMETER ::   jpr_cal    = 29            ! calving
95   INTEGER, PARAMETER ::   jpr_taum   = 30            ! wind stress module
96   INTEGER, PARAMETER ::   jpr_co2    = 31
97   INTEGER, PARAMETER ::   jpr_topm   = 32            ! topmeltn
98   INTEGER, PARAMETER ::   jpr_botm   = 33            ! botmeltn
99   INTEGER, PARAMETER ::   jprcv      = 33            ! total number of fields received
100
101   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction
102   INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
103   INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
104   INTEGER, PARAMETER ::   jps_tmix   =  4            ! mixed temperature (ocean+ice)
105   INTEGER, PARAMETER ::   jps_albice =  5            ! ice   albedo
106   INTEGER, PARAMETER ::   jps_albmix =  6            ! mixed albedo
107   INTEGER, PARAMETER ::   jps_hice   =  7            ! ice  thickness
108   INTEGER, PARAMETER ::   jps_hsnw   =  8            ! snow thickness
109   INTEGER, PARAMETER ::   jps_ocx1   =  9            ! ocean current on grid 1
110   INTEGER, PARAMETER ::   jps_ocy1   = 10            !
111   INTEGER, PARAMETER ::   jps_ocz1   = 11            !
112   INTEGER, PARAMETER ::   jps_ivx1   = 12            ! ice   current on grid 1
113   INTEGER, PARAMETER ::   jps_ivy1   = 13            !
114   INTEGER, PARAMETER ::   jps_ivz1   = 14            !
115   INTEGER, PARAMETER ::   jps_co2    = 15
116   INTEGER, PARAMETER ::   jpsnd      = 15            ! total number of fields sended
117
118   !                                                         !!** namelist namsbc_cpl **
119   TYPE ::   FLD_C
120      CHARACTER(len = 32) ::   cldes                  ! desciption of the coupling strategy
121      CHARACTER(len = 32) ::   clcat                  ! multiple ice categories strategy
122      CHARACTER(len = 32) ::   clvref                 ! reference of vector ('spherical' or 'cartesian')
123      CHARACTER(len = 32) ::   clvor                  ! orientation of vector fields ('eastward-northward' or 'local grid')
124      CHARACTER(len = 32) ::   clvgrd                 ! grids on which is located the vector fields
125   END TYPE FLD_C
126   ! Send to the atmosphere                           !
127   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                       
128   ! Received from the atmosphere                     !
129   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
130   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                       
131
132   TYPE ::   DYNARR     
133      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
134   END TYPE DYNARR
135
136   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                      ! all fields recieved from the atmosphere
137
138   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
139
140   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
141
142#if ! defined key_lim2   &&   ! defined key_lim3
143   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice,fr1_i0,fr2_i0          ! jpi, jpj
144   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl)
145#endif
146
147#if defined key_cice
148   INTEGER, PARAMETER ::   jpl = ncat
149#elif ! defined key_lim2   &&   ! defined key_lim3
150   INTEGER, PARAMETER ::   jpl = 1 
151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice
152   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice
153#endif
154
155#if ! defined key_lim3   &&  ! defined key_cice
156   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i
157#endif
158
159#if ! defined key_lim3
160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  ht_i, ht_s
161#endif
162
163#if ! defined key_cice
164   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  topmelt, botmelt
165#endif
166
167   !! Substitution
168#  include "vectopt_loop_substitute.h90"
169   !!----------------------------------------------------------------------
170   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
171   !! $Id$
172   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
173   !!----------------------------------------------------------------------
174
175CONTAINS
176 
177   INTEGER FUNCTION sbc_cpl_alloc()
178      !!----------------------------------------------------------------------
179      !!             ***  FUNCTION sbc_cpl_alloc  ***
180      !!----------------------------------------------------------------------
181      INTEGER :: ierr(4),jn
182      !!----------------------------------------------------------------------
183      ierr(:) = 0
184      !
185      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
186      !
187#if ! defined key_lim2 && ! defined key_lim3
188      ! quick patch to be able to run the coupled model without sea-ice...
189      ALLOCATE( u_ice(jpi,jpj) , fr1_i0(jpi,jpj) , tn_ice (jpi,jpj,1) ,     &
190                v_ice(jpi,jpj) , fr2_i0(jpi,jpj) , alb_ice(jpi,jpj,1),      &
191                emp_ice(jpi,jpj) , qns_ice(jpi,jpj,1) , dqns_ice(jpi,jpj,1) , STAT=ierr(2) )
192#endif
193
194#if ! defined key_lim3 && ! defined key_cice
195      ALLOCATE( a_i(jpi,jpj,jpl) , STAT=ierr(3) )
196#endif
197
198#if defined key_cice || defined key_lim2
199      ALLOCATE( ht_i(jpi,jpj,jpl) , ht_s(jpi,jpj,jpl) , STAT=ierr(4) )
200#endif
201      sbc_cpl_alloc = MAXVAL( ierr )
202      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
203      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
204      !
205   END FUNCTION sbc_cpl_alloc
206
207
208   SUBROUTINE sbc_cpl_init( k_ice )     
209      !!----------------------------------------------------------------------
210      !!             ***  ROUTINE sbc_cpl_init  ***
211      !!
212      !! ** Purpose :   Initialisation of send and recieved information from
213      !!                the atmospheric component
214      !!
215      !! ** Method  : * Read namsbc_cpl namelist
216      !!              * define the receive interface
217      !!              * define the send    interface
218      !!              * initialise the OASIS coupler
219      !!----------------------------------------------------------------------
220      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
221      !!
222      INTEGER ::   jn   ! dummy loop index
223      INTEGER ::   ios  ! Local integer output status for namelist read
224      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
225      !!
226      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,   &
227         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,   &
228         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx  , sn_rcv_co2
229      !!---------------------------------------------------------------------
230      !
231      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
232      !
233      CALL wrk_alloc( jpi,jpj, zacs, zaos )
234
235      ! ================================ !
236      !      Namelist informations       !
237      ! ================================ !
238
239      REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
240      READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
241901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
242
243      REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
244      READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
245902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
246      WRITE ( numond, namsbc_cpl )
247
248      IF(lwp) THEN                        ! control print
249         WRITE(numout,*)
250         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
251         WRITE(numout,*)'~~~~~~~~~~~~'
252         WRITE(numout,*)'  received fields (mutiple ice categogies)'
253         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
254         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
255         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
256         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
257         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
258         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
259         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
260         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
261         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
262         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
263         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
264         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
265         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
266         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
267         WRITE(numout,*)'  sent fields (multiple ice categories)'
268         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
269         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
270         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
271         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
272         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
273         WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
274         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
275         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
276      ENDIF
277
278      !                                   ! allocate sbccpl arrays
279      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
280     
281      ! ================================ !
282      !   Define the receive interface   !
283      ! ================================ !
284      nrcvinfo(:) = OASIS_idle   ! needed by nrcvinfo(jpr_otx1) if we do not receive ocean stress
285
286      ! for each field: define the OASIS name                              (srcv(:)%clname)
287      !                 define receive or not from the namelist parameters (srcv(:)%laction)
288      !                 define the north fold type of lbc                  (srcv(:)%nsgn)
289
290      ! default definitions of srcv
291      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
292
293      !                                                      ! ------------------------- !
294      !                                                      ! ice and ocean wind stress !   
295      !                                                      ! ------------------------- !
296      !                                                           ! Name
297      srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
298      srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
299      srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
300      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
301      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
302      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
303      !
304      srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
305      srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
306      srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
307      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
308      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
309      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
310      !
311      ! Vectors: change of sign at north fold ONLY if on the local grid
312      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
313     
314      !                                                           ! Set grid and action
315      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
316      CASE( 'T' ) 
317         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
318         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
319         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
320      CASE( 'U,V' ) 
321         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
322         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
323         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
324         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
325         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
326      CASE( 'U,V,T' )
327         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
328         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
329         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
330         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
331         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
332      CASE( 'U,V,I' )
333         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
334         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
335         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
336         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
337         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
338      CASE( 'U,V,F' )
339         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
340         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
341         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
342         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
343         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
344      CASE( 'T,I' ) 
345         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
346         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
347         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
348         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
349      CASE( 'T,F' ) 
350         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
351         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
352         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
353         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
354      CASE( 'T,U,V' )
355         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
356         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
357         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
358         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
359         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
360      CASE default   
361         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
362      END SELECT
363      !
364      IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
365         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE. 
366      !
367      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
368            srcv(jpr_otx2:jpr_otz2)%laction = .FALSE. 
369            srcv(jpr_itx2:jpr_itz2)%laction = .FALSE. 
370            srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid   ! not needed but cleaner...
371            srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid   ! not needed but cleaner...
372      ENDIF
373      !
374      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
375         srcv(jpr_itz1:jpr_itz2)%laction = .FALSE.    ! ice components not received (itx1 and ity1 used later)
376         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
377         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
378      ENDIF
379       
380      !                                                      ! ------------------------- !
381      !                                                      !    freshwater budget      !   E-P
382      !                                                      ! ------------------------- !
383      ! we suppose that atmosphere modele do not make the difference between precipiration (liquide or solid)
384      ! over ice of free ocean within the same atmospheric cell.cd
385      srcv(jpr_rain)%clname = 'OTotRain'      ! Rain = liquid precipitation
386      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
387      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
388      srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation
389      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation
390      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
391      srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
392      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
393      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
394      CASE( 'conservative'  )   ;   srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
395      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
396      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
397      END SELECT
398
399      !                                                      ! ------------------------- !
400      !                                                      !     Runoffs & Calving     !   
401      !                                                      ! ------------------------- !
402      srcv(jpr_rnf   )%clname = 'O_Runoff'   ;   IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' )   srcv(jpr_rnf)%laction = .TRUE.
403! This isn't right - really just want ln_rnf_emp changed
404!                                                 IF( TRIM( sn_rcv_rnf%cldes ) == 'climato' )   THEN   ;   ln_rnf = .TRUE.
405!                                                 ELSE                                                 ;   ln_rnf = .FALSE.
406!                                                 ENDIF
407      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
408
409      !                                                      ! ------------------------- !
410      !                                                      !    non solar radiation    !   Qns
411      !                                                      ! ------------------------- !
412      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
413      srcv(jpr_qnsice)%clname = 'O_QnsIce'
414      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
415      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
416      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
417      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
418      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
419      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
420      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
421      END SELECT
422      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
423         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
424      !                                                      ! ------------------------- !
425      !                                                      !    solar radiation        !   Qsr
426      !                                                      ! ------------------------- !
427      srcv(jpr_qsroce)%clname = 'O_QsrOce'
428      srcv(jpr_qsrice)%clname = 'O_QsrIce'
429      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
430      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
431      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
432      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
433      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
434      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
435      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
436      END SELECT
437      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
438         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
439      !                                                      ! ------------------------- !
440      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
441      !                                                      ! ------------------------- !
442      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
443      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
444      !
445      ! non solar sensitivity mandatory for LIM ice model
446      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4) &
447         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
448      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
449      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
450         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
451      !                                                      ! ------------------------- !
452      !                                                      !    Ice Qsr penetration    !   
453      !                                                      ! ------------------------- !
454      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer
455      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
456      ! Coupled case: since cloud cover is not received from atmosphere
457      !               ===> defined as constant value -> definition done in sbc_cpl_init
458      fr1_i0(:,:) = 0.18
459      fr2_i0(:,:) = 0.82
460      !                                                      ! ------------------------- !
461      !                                                      !      10m wind module      !   
462      !                                                      ! ------------------------- !
463      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
464      !
465      !                                                      ! ------------------------- !
466      !                                                      !   wind stress module      !   
467      !                                                      ! ------------------------- !
468      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
469      lhftau = srcv(jpr_taum)%laction
470
471      !                                                      ! ------------------------- !
472      !                                                      !      Atmospheric CO2      !
473      !                                                      ! ------------------------- !
474      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
475      !                                                      ! ------------------------- !
476      !                                                      !   topmelt and botmelt     !   
477      !                                                      ! ------------------------- !
478      srcv(jpr_topm )%clname = 'OTopMlt'
479      srcv(jpr_botm )%clname = 'OBotMlt'
480      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
481         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
482            srcv(jpr_topm:jpr_botm)%nct = jpl
483         ELSE
484            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
485         ENDIF
486         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
487      ENDIF
488
489      ! Allocate all parts of frcv used for received fields
490      DO jn = 1, jprcv
491         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
492      END DO
493      ! Allocate taum part of frcv which is used even when not received as coupling field
494      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jn)%nct) )
495
496      ! ================================ !
497      !     Define the send interface    !
498      ! ================================ !
499      ! for each field: define the OASIS name                           (ssnd(:)%clname)
500      !                 define send or not from the namelist parameters (ssnd(:)%laction)
501      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
502     
503      ! default definitions of nsnd
504      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
505         
506      !                                                      ! ------------------------- !
507      !                                                      !    Surface temperature    !
508      !                                                      ! ------------------------- !
509      ssnd(jps_toce)%clname = 'O_SSTSST'
510      ssnd(jps_tice)%clname = 'O_TepIce'
511      ssnd(jps_tmix)%clname = 'O_TepMix'
512      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
513      CASE( 'none'         )       ! nothing to do
514      CASE( 'oce only'             )   ;   ssnd(   jps_toce             )%laction = .TRUE.
515      CASE( 'weighted oce and ice' )
516         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
517         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
518      CASE( 'mixed oce-ice'        )   ;   ssnd(   jps_tmix             )%laction = .TRUE.
519      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
520      END SELECT
521     
522      !                                                      ! ------------------------- !
523      !                                                      !          Albedo           !
524      !                                                      ! ------------------------- !
525      ssnd(jps_albice)%clname = 'O_AlbIce' 
526      ssnd(jps_albmix)%clname = 'O_AlbMix'
527      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
528      CASE( 'none'          )       ! nothing to do
529      CASE( 'weighted ice'  )   ;   ssnd(jps_albice)%laction = .TRUE.
530      CASE( 'mixed oce-ice' )   ;   ssnd(jps_albmix)%laction = .TRUE.
531      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
532      END SELECT
533      !
534      ! Need to calculate oceanic albedo if
535      !     1. sending mixed oce-ice albedo or
536      !     2. receiving mixed oce-ice solar radiation
537      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
538         CALL albedo_oce( zaos, zacs )
539         ! Due to lack of information on nebulosity : mean clear/overcast sky
540         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
541      ENDIF
542
543      !                                                      ! ------------------------- !
544      !                                                      !  Ice fraction & Thickness !
545      !                                                      ! ------------------------- !
546      ssnd(jps_fice)%clname = 'OIceFrc'
547      ssnd(jps_hice)%clname = 'OIceTck'
548      ssnd(jps_hsnw)%clname = 'OSnwTck'
549      IF( k_ice /= 0 ) THEN
550         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
551! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
552         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
553      ENDIF
554
555      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
556      CASE( 'none'         )       ! nothing to do
557      CASE( 'ice and snow' ) 
558         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
559         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
560            ssnd(jps_hice:jps_hsnw)%nct = jpl
561         ELSE
562            IF ( jpl > 1 ) THEN
563CALL ctl_stop( 'sbc_cpl_init: use weighted ice and snow option for sn_snd_thick%cldes if not exchanging category fields' )
564            ENDIF
565         ENDIF
566      CASE ( 'weighted ice and snow' ) 
567         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
568         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
569      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
570      END SELECT
571
572      !                                                      ! ------------------------- !
573      !                                                      !      Surface current      !
574      !                                                      ! ------------------------- !
575      !        ocean currents              !            ice velocities
576      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
577      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
578      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
579      !
580      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
581
582      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
583         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
584      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
585         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
586         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
587      ENDIF
588      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
589      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
590      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
591      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
592      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
593      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
594      CASE( 'weighted oce and ice' )   !   nothing to do
595      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
596      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
597      END SELECT
598
599      !                                                      ! ------------------------- !
600      !                                                      !          CO2 flux         !
601      !                                                      ! ------------------------- !
602      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
603      !
604      ! ================================ !
605      !   initialisation of the coupler  !
606      ! ================================ !
607
608      CALL cpl_prism_define(jprcv, jpsnd)           
609      !
610      IF( ln_dm2dc .AND. ( cpl_prism_freq( jpr_qsroce ) + cpl_prism_freq( jpr_qsrmix ) /= 86400 ) )   &
611         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
612
613      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
614      !
615      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
616      !
617   END SUBROUTINE sbc_cpl_init
618
619
620   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
621      !!----------------------------------------------------------------------
622      !!             ***  ROUTINE sbc_cpl_rcv  ***
623      !!
624      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
625      !!                provide the ocean heat and freshwater fluxes.
626      !!
627      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
628      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
629      !!                to know if the field was really received or not
630      !!
631      !!              --> If ocean stress was really received:
632      !!
633      !!                  - transform the received ocean stress vector from the received
634      !!                 referential and grid into an atmosphere-ocean stress in
635      !!                 the (i,j) ocean referencial and at the ocean velocity point.
636      !!                    The received stress are :
637      !!                     - defined by 3 components (if cartesian coordinate)
638      !!                            or by 2 components (if spherical)
639      !!                     - oriented along geographical   coordinate (if eastward-northward)
640      !!                            or  along the local grid coordinate (if local grid)
641      !!                     - given at U- and V-point, resp.   if received on 2 grids
642      !!                            or at T-point               if received on 1 grid
643      !!                    Therefore and if necessary, they are successively
644      !!                  processed in order to obtain them
645      !!                     first  as  2 components on the sphere
646      !!                     second as  2 components oriented along the local grid
647      !!                     third  as  2 components on the U,V grid
648      !!
649      !!              -->
650      !!
651      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
652      !!             and total ocean freshwater fluxes 
653      !!
654      !! ** Method  :   receive all fields from the atmosphere and transform
655      !!              them into ocean surface boundary condition fields
656      !!
657      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
658      !!                        taum, wndm   wind stres and wind speed module at T-point
659      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
660      !!                                     and the latent heat flux of solid precip. melting
661      !!                        qsr          solar ocean heat fluxes   (ocean only case)
662      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
663      !!----------------------------------------------------------------------
664      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
665      INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation
666      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
667      !!
668      LOGICAL ::    llnewtx, llnewtau      ! update wind stress components and module??
669      INTEGER  ::   ji, jj, jn             ! dummy loop indices
670      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
671      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
672      REAL(wp) ::   zcoef                  ! temporary scalar
673      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
674      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
675      REAL(wp) ::   zzx, zzy               ! temporary variables
676      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
677      !!----------------------------------------------------------------------
678      !
679      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
680      !
681      CALL wrk_alloc( jpi,jpj, ztx, zty )
682
683      IF( kt == nit000 )   CALL sbc_cpl_init( k_ice )          ! initialisation
684
685      !                                                 ! Receive all the atmos. fields (including ice information)
686      isec = ( kt - nit000 ) * NINT( rdttra(1) )             ! date of exchanges
687      DO jn = 1, jprcv                                       ! received fields sent by the atmosphere
688         IF( srcv(jn)%laction )   CALL cpl_prism_rcv( jn, isec, frcv(jn)%z3, nrcvinfo(jn) )
689      END DO
690
691      !                                                      ! ========================= !
692      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
693         !                                                   ! ========================= !
694         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
695         ! => need to be done only when we receive the field
696         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
697            !
698            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
699               !                                                       ! (cartesian to spherical -> 3 to 2 components)
700               !
701               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
702                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
703               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
704               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
705               !
706               IF( srcv(jpr_otx2)%laction ) THEN
707                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
708                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
709                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
710                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
711               ENDIF
712               !
713            ENDIF
714            !
715            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
716               !                                                       ! (geographical to local grid -> rotate the components)
717               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
718               IF( srcv(jpr_otx2)%laction ) THEN
719                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
720               ELSE 
721                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
722               ENDIF
723               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
724               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
725            ENDIF
726            !                             
727            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
728               DO jj = 2, jpjm1                                          ! T ==> (U,V)
729                  DO ji = fs_2, fs_jpim1   ! vector opt.
730                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
731                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
732                  END DO
733               END DO
734               CALL lbc_lnk( frcv(jpr_otx1)%z3(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
735            ENDIF
736            llnewtx = .TRUE.
737         ELSE
738            llnewtx = .FALSE.
739         ENDIF
740         !                                                   ! ========================= !
741      ELSE                                                   !   No dynamical coupling   !
742         !                                                   ! ========================= !
743         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
744         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
745         llnewtx = .TRUE.
746         !
747      ENDIF
748     
749      !                                                      ! ========================= !
750      !                                                      !    wind stress module     !   (taum)
751      !                                                      ! ========================= !
752      !
753      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
754         ! => need to be done only when otx1 was changed
755         IF( llnewtx ) THEN
756!CDIR NOVERRCHK
757            DO jj = 2, jpjm1
758!CDIR NOVERRCHK
759               DO ji = fs_2, fs_jpim1   ! vect. opt.
760                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
761                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
762                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
763               END DO
764            END DO
765            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
766            llnewtau = .TRUE.
767         ELSE
768            llnewtau = .FALSE.
769         ENDIF
770      ELSE
771         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
772         ! Stress module can be negative when received (interpolation problem)
773         IF( llnewtau ) THEN
774            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
775         ENDIF
776      ENDIF
777     
778      !                                                      ! ========================= !
779      !                                                      !      10 m wind speed      !   (wndm)
780      !                                                      ! ========================= !
781      !
782      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
783         ! => need to be done only when taumod was changed
784         IF( llnewtau ) THEN
785            zcoef = 1. / ( zrhoa * zcdrag ) 
786!CDIR NOVERRCHK
787            DO jj = 1, jpj
788!CDIR NOVERRCHK
789               DO ji = 1, jpi 
790                  wndm(ji,jj) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
791               END DO
792            END DO
793         ENDIF
794      ELSE
795         IF ( nrcvinfo(jpr_w10m) == OASIS_Rcv ) wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
796      ENDIF
797
798      ! u(v)tau and taum will be modified by ice model
799      ! -> need to be reset before each call of the ice/fsbc     
800      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
801         !
802         utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
803         vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
804         taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
805         CALL iom_put( "taum_oce", taum )   ! output wind stress module
806         
807      ENDIF
808
809#if defined key_cpl_carbon_cycle
810      !                                                              ! atmosph. CO2 (ppm)
811      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
812#endif
813
814      !                                                      ! ========================= !
815      IF( k_ice <= 1 ) THEN                                  !  heat & freshwater fluxes ! (Ocean only case)
816         !                                                   ! ========================= !
817         !
818         !                                                       ! total freshwater fluxes over the ocean (emp)
819         SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
820         CASE( 'conservative' )
821            emp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
822         CASE( 'oce only', 'oce and ice' )
823            emp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
824         CASE default
825            CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
826         END SELECT
827         !
828         !                                                        ! runoffs and calving (added in emp)
829         IF( srcv(jpr_rnf)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
830         IF( srcv(jpr_cal)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_cal)%z3(:,:,1)
831         !
832!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
833!!gm                                       at least should be optional...
834!!         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN     ! add to the total freshwater budget
835!!            ! remove negative runoff
836!!            zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
837!!            zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
838!!            IF( lk_mpp )   CALL mpp_sum( zcumulpos )   ! sum over the global domain
839!!            IF( lk_mpp )   CALL mpp_sum( zcumulneg )
840!!            IF( zcumulpos /= 0. ) THEN                 ! distribute negative runoff on positive runoff grid points
841!!               zcumulneg = 1.e0 + zcumulneg / zcumulpos
842!!               frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
843!!            ENDIF     
844!!            ! add runoff to e-p
845!!            emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
846!!         ENDIF
847!!gm  end of internal cooking
848         !
849         !                                                       ! non solar heat flux over the ocean (qns)
850         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
851         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
852         ! add the latent heat of solid precip. melting
853         IF( srcv(jpr_snow  )%laction )   THEN                         ! update qns over the free ocean with:
854              qns(:,:) = qns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus  & ! energy for melting solid precipitation over the free ocean
855           &           - emp(:,:) * sst_m(:,:) * rcp                   ! remove heat content due to mass flux (assumed to be at SST)
856         ENDIF
857
858         !                                                       ! solar flux over the ocean          (qsr)
859         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
860         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
861         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
862         !
863 
864      ENDIF
865      !
866      CALL wrk_dealloc( jpi,jpj, ztx, zty )
867      !
868      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
869      !
870   END SUBROUTINE sbc_cpl_rcv
871   
872
873   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
874      !!----------------------------------------------------------------------
875      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
876      !!
877      !! ** Purpose :   provide the stress over sea-ice in coupled mode
878      !!
879      !! ** Method  :   transform the received stress from the atmosphere into
880      !!             an atmosphere-ice stress in the (i,j) ocean referencial
881      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
882      !!                'C'-grid : i- (j-) components given at U- (V-) point
883      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
884      !!
885      !!                The received stress are :
886      !!                 - defined by 3 components (if cartesian coordinate)
887      !!                        or by 2 components (if spherical)
888      !!                 - oriented along geographical   coordinate (if eastward-northward)
889      !!                        or  along the local grid coordinate (if local grid)
890      !!                 - given at U- and V-point, resp.   if received on 2 grids
891      !!                        or at a same point (T or I) if received on 1 grid
892      !!                Therefore and if necessary, they are successively
893      !!             processed in order to obtain them
894      !!                 first  as  2 components on the sphere
895      !!                 second as  2 components oriented along the local grid
896      !!                 third  as  2 components on the cp_ice_msh point
897      !!
898      !!                Except in 'oce and ice' case, only one vector stress field
899      !!             is received. It has already been processed in sbc_cpl_rcv
900      !!             so that it is now defined as (i,j) components given at U-
901      !!             and V-points, respectively. Therefore, only the third
902      !!             transformation is done and only if the ice-grid is a 'I'-grid.
903      !!
904      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
905      !!----------------------------------------------------------------------
906      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
907      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
908      !!
909      INTEGER ::   ji, jj                          ! dummy loop indices
910      INTEGER ::   itx                             ! index of taux over ice
911      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
912      !!----------------------------------------------------------------------
913      !
914      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
915      !
916      CALL wrk_alloc( jpi,jpj, ztx, zty )
917
918      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
919      ELSE                                ;   itx =  jpr_otx1
920      ENDIF
921
922      ! do something only if we just received the stress from atmosphere
923      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
924
925         !                                                      ! ======================= !
926         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
927            !                                                   ! ======================= !
928           
929            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
930               !                                                       ! (cartesian to spherical -> 3 to 2 components)
931               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
932                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
933               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
934               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
935               !
936               IF( srcv(jpr_itx2)%laction ) THEN
937                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
938                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
939                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
940                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
941               ENDIF
942               !
943            ENDIF
944            !
945            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
946               !                                                       ! (geographical to local grid -> rotate the components)
947               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
948               IF( srcv(jpr_itx2)%laction ) THEN
949                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
950               ELSE
951                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
952               ENDIF
953               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
954               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
955            ENDIF
956            !                                                   ! ======================= !
957         ELSE                                                   !     use ocean stress    !
958            !                                                   ! ======================= !
959            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
960            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
961            !
962         ENDIF
963
964         !                                                      ! ======================= !
965         !                                                      !     put on ice grid     !
966         !                                                      ! ======================= !
967         !   
968         !                                                  j+1   j     -----V---F
969         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
970         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
971         !                                                               |       |
972         !                                                   j    j-1   -I-------|
973         !                                               (for I)         |       |
974         !                                                              i-1  i   i
975         !                                                               i      i+1 (for I)
976         SELECT CASE ( cp_ice_msh )
977            !
978         CASE( 'I' )                                         ! B-grid ==> I
979            SELECT CASE ( srcv(jpr_itx1)%clgrid )
980            CASE( 'U' )
981               DO jj = 2, jpjm1                                   ! (U,V) ==> I
982                  DO ji = 2, jpim1   ! NO vector opt.
983                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
984                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
985                  END DO
986               END DO
987            CASE( 'F' )
988               DO jj = 2, jpjm1                                   ! F ==> I
989                  DO ji = 2, jpim1   ! NO vector opt.
990                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
991                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
992                  END DO
993               END DO
994            CASE( 'T' )
995               DO jj = 2, jpjm1                                   ! T ==> I
996                  DO ji = 2, jpim1   ! NO vector opt.
997                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
998                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
999                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
1000                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1001                  END DO
1002               END DO
1003            CASE( 'I' )
1004               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
1005               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1006            END SELECT
1007            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1008               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1009            ENDIF
1010            !
1011         CASE( 'F' )                                         ! B-grid ==> F
1012            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1013            CASE( 'U' )
1014               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1015                  DO ji = fs_2, fs_jpim1   ! vector opt.
1016                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1017                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1018                  END DO
1019               END DO
1020            CASE( 'I' )
1021               DO jj = 2, jpjm1                                   ! I ==> F
1022                  DO ji = 2, jpim1   ! NO vector opt.
1023                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1024                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1025                  END DO
1026               END DO
1027            CASE( 'T' )
1028               DO jj = 2, jpjm1                                   ! T ==> F
1029                  DO ji = 2, jpim1   ! NO vector opt.
1030                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1031                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1032                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1033                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1034                  END DO
1035               END DO
1036            CASE( 'F' )
1037               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1038               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1039            END SELECT
1040            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1041               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1042            ENDIF
1043            !
1044         CASE( 'C' )                                         ! C-grid ==> U,V
1045            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1046            CASE( 'U' )
1047               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1048               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1049            CASE( 'F' )
1050               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1051                  DO ji = fs_2, fs_jpim1   ! vector opt.
1052                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1053                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1054                  END DO
1055               END DO
1056            CASE( 'T' )
1057               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1058                  DO ji = fs_2, fs_jpim1   ! vector opt.
1059                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1060                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1061                  END DO
1062               END DO
1063            CASE( 'I' )
1064               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1065                  DO ji = 2, jpim1   ! NO vector opt.
1066                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1067                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1068                  END DO
1069               END DO
1070            END SELECT
1071            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1072               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1073            ENDIF
1074         END SELECT
1075
1076      ENDIF
1077      !   
1078      CALL wrk_dealloc( jpi,jpj, ztx, zty )
1079      !
1080      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
1081      !
1082   END SUBROUTINE sbc_cpl_ice_tau
1083   
1084
1085   SUBROUTINE sbc_cpl_ice_flx( p_frld  , palbi   , psst    , pist    )
1086      !!----------------------------------------------------------------------
1087      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1088      !!
1089      !! ** Purpose :   provide the heat and freshwater fluxes of the
1090      !!              ocean-ice system.
1091      !!
1092      !! ** Method  :   transform the fields received from the atmosphere into
1093      !!             surface heat and fresh water boundary condition for the
1094      !!             ice-ocean system. The following fields are provided:
1095      !!              * total non solar, solar and freshwater fluxes (qns_tot,
1096      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1097      !!             NB: emp_tot include runoffs and calving.
1098      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1099      !!             emp_ice = sublimation - solid precipitation as liquid
1100      !!             precipitation are re-routed directly to the ocean and
1101      !!             runoffs and calving directly enter the ocean.
1102      !!              * solid precipitation (sprecip), used to add to qns_tot
1103      !!             the heat lost associated to melting solid precipitation
1104      !!             over the ocean fraction.
1105      !!       ===>> CAUTION here this changes the net heat flux received from
1106      !!             the atmosphere
1107      !!
1108      !!                  - the fluxes have been separated from the stress as
1109      !!                 (a) they are updated at each ice time step compare to
1110      !!                 an update at each coupled time step for the stress, and
1111      !!                 (b) the conservative computation of the fluxes over the
1112      !!                 sea-ice area requires the knowledge of the ice fraction
1113      !!                 after the ice advection and before the ice thermodynamics,
1114      !!                 so that the stress is updated before the ice dynamics
1115      !!                 while the fluxes are updated after it.
1116      !!
1117      !! ** Action  :   update at each nf_ice time step:
1118      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1119      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1120      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
1121      !!                   emp_ice            ice sublimation - solid precipitation over the ice
1122      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
1123      !!                   sprecip             solid precipitation over the ocean 
1124      !!----------------------------------------------------------------------
1125      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1126      ! optional arguments, used only in 'mixed oce-ice' case
1127      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1128      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celcius]
1129      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature     [Kelvin]
1130      !
1131      INTEGER ::   jl   ! dummy loop index
1132      REAL(wp), POINTER, DIMENSION(:,:) ::   zcptn, ztmp, zicefr
1133      !!----------------------------------------------------------------------
1134      !
1135      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1136      !
1137      CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr )
1138
1139      zicefr(:,:) = 1.- p_frld(:,:)
1140      zcptn(:,:) = rcp * sst_m(:,:)
1141      !
1142      !                                                      ! ========================= !
1143      !                                                      !    freshwater budget      !   (emp)
1144      !                                                      ! ========================= !
1145      !
1146      !                                                           ! total Precipitations - total Evaporation (emp_tot)
1147      !                                                           ! solid precipitation  - sublimation       (emp_ice)
1148      !                                                           ! solid Precipitation                      (sprecip)
1149      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1150      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1151         sprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                 ! May need to ensure positive here
1152         tprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + sprecip (:,:) ! May need to ensure positive here
1153         emp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - tprecip(:,:)
1154         emp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1155                           CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation
1156         IF( lk_diaar5 )   CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
1157         ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1158                           CALL iom_put( 'evap_ao_cea'  , ztmp                            )   ! ice-free oce evap (cell average)
1159         IF( lk_diaar5 )   CALL iom_put( 'hflx_evap_cea', ztmp(:,:         ) * zcptn(:,:) )   ! heat flux from from evap (cell ave)
1160      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1161         emp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1162         emp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
1163         sprecip(:,:) = - frcv(jpr_semp)%z3(:,:,1) + frcv(jpr_ievp)%z3(:,:,1)
1164      END SELECT
1165
1166      CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
1167      CALL iom_put( 'snow_ao_cea', sprecip(:,:         ) * p_frld(:,:)    )   ! Snow        over ice-free ocean  (cell average)
1168      CALL iom_put( 'snow_ai_cea', sprecip(:,:         ) * zicefr(:,:)    )   ! Snow        over sea-ice         (cell average)
1169      CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
1170      !   
1171      !                                                           ! runoffs and calving (put in emp_tot)
1172      IF( srcv(jpr_rnf)%laction ) THEN
1173         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)
1174                           CALL iom_put( 'runoffs'      , frcv(jpr_rnf)%z3(:,:,1)              )   ! rivers
1175         IF( lk_diaar5 )   CALL iom_put( 'hflx_rnf_cea' , frcv(jpr_rnf)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from rivers
1176      ENDIF
1177      IF( srcv(jpr_cal)%laction ) THEN
1178         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1179         CALL iom_put( 'calving', frcv(jpr_cal)%z3(:,:,1) )
1180      ENDIF
1181      !
1182!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
1183!!gm                                       at least should be optional...
1184!!       ! remove negative runoff                            ! sum over the global domain
1185!!       zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
1186!!       zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
1187!!       IF( lk_mpp )   CALL mpp_sum( zcumulpos )
1188!!       IF( lk_mpp )   CALL mpp_sum( zcumulneg )
1189!!       IF( zcumulpos /= 0. ) THEN                          ! distribute negative runoff on positive runoff grid points
1190!!          zcumulneg = 1.e0 + zcumulneg / zcumulpos
1191!!          frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
1192!!       ENDIF     
1193!!       emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)   ! add runoff to e-p
1194!!
1195!!gm  end of internal cooking
1196
1197      !                                                      ! ========================= !
1198      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1199      !                                                      ! ========================= !
1200      CASE( 'oce only' )                                     ! the required field is directly provided
1201         qns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
1202      CASE( 'conservative' )                                      ! the required fields are directly provided
1203         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1204         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1205            qns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1206         ELSE
1207            ! Set all category values equal for the moment
1208            DO jl=1,jpl
1209               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1210            ENDDO
1211         ENDIF
1212      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
1213         qns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1214         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1215            DO jl=1,jpl
1216               qns_tot(:,:   ) = qns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1217               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1218            ENDDO
1219         ELSE
1220            DO jl=1,jpl
1221               qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1222               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1223            ENDDO
1224         ENDIF
1225      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
1226! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1227         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1228         qns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1229            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1230            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
1231      END SELECT
1232      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus
1233      qns_tot(:,:) = qns_tot(:,:)                         &            ! qns_tot update over free ocean with:
1234         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
1235         &          - (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
1236         &             - emp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
1237      IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1238!!gm
1239!!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in
1240!!    the flux that enter the ocean....
1241!!    moreover 1 - it is not diagnose anywhere....
1242!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1243!!
1244!! similar job should be done for snow and precipitation temperature
1245      !                                     
1246      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting
1247         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting
1248         qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)
1249         IF( lk_diaar5 )   CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1250      ENDIF
1251
1252      !                                                      ! ========================= !
1253      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1254      !                                                      ! ========================= !
1255      CASE( 'oce only' )
1256         qsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
1257      CASE( 'conservative' )
1258         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1259         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1260            qsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1261         ELSE
1262            ! Set all category values equal for the moment
1263            DO jl=1,jpl
1264               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1265            ENDDO
1266         ENDIF
1267         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1268         qsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1269      CASE( 'oce and ice' )
1270         qsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1271         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1272            DO jl=1,jpl
1273               qsr_tot(:,:   ) = qsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1274               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1275            ENDDO
1276         ELSE
1277            DO jl=1,jpl
1278               qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1279               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1280            ENDDO
1281         ENDIF
1282      CASE( 'mixed oce-ice' )
1283         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1284! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1285!       Create solar heat flux over ice using incoming solar heat flux and albedos
1286!       ( see OASIS3 user guide, 5th edition, p39 )
1287         qsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1288            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1289            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1290      END SELECT
1291      IF( ln_dm2dc ) THEN   ! modify qsr to include the diurnal cycle
1292         qsr_tot(:,:  ) = sbc_dcy( qsr_tot(:,:  ) )
1293         DO jl=1,jpl
1294            qsr_ice(:,:,jl) = sbc_dcy( qsr_ice(:,:,jl) )
1295         ENDDO
1296      ENDIF
1297
1298      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )
1299      CASE ('coupled')
1300         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1301            dqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
1302         ELSE
1303            ! Set all category values equal for the moment
1304            DO jl=1,jpl
1305               dqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
1306            ENDDO
1307         ENDIF
1308      END SELECT
1309
1310      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )
1311      CASE ('coupled')
1312         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
1313         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
1314      END SELECT
1315
1316      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
1317      !
1318      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
1319      !
1320   END SUBROUTINE sbc_cpl_ice_flx
1321   
1322   
1323   SUBROUTINE sbc_cpl_snd( kt )
1324      !!----------------------------------------------------------------------
1325      !!             ***  ROUTINE sbc_cpl_snd  ***
1326      !!
1327      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
1328      !!
1329      !! ** Method  :   send to the atmosphere through a call to cpl_prism_snd
1330      !!              all the needed fields (as defined in sbc_cpl_init)
1331      !!----------------------------------------------------------------------
1332      INTEGER, INTENT(in) ::   kt
1333      !
1334      INTEGER ::   ji, jj, jl   ! dummy loop indices
1335      INTEGER ::   isec, info   ! local integer
1336      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
1337      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
1338      !!----------------------------------------------------------------------
1339      !
1340      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
1341      !
1342      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1343      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1344
1345      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
1346
1347      zfr_l(:,:) = 1.- fr_i(:,:)
1348
1349      !                                                      ! ------------------------- !
1350      !                                                      !    Surface temperature    !   in Kelvin
1351      !                                                      ! ------------------------- !
1352      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
1353         SELECT CASE( sn_snd_temp%cldes)
1354         CASE( 'oce only'             )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0
1355         CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:)   
1356            SELECT CASE( sn_snd_temp%clcat )
1357            CASE( 'yes' )   
1358               ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1359            CASE( 'no' )
1360               ztmp3(:,:,:) = 0.0
1361               DO jl=1,jpl
1362                  ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
1363               ENDDO
1364            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1365            END SELECT
1366         CASE( 'mixed oce-ice'        )   
1367            ztmp1(:,:) = ( tsn(:,:,1,1) + rt0 ) * zfr_l(:,:) 
1368            DO jl=1,jpl
1369               ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
1370            ENDDO
1371         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
1372         END SELECT
1373         IF( ssnd(jps_toce)%laction )   CALL cpl_prism_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1374         IF( ssnd(jps_tice)%laction )   CALL cpl_prism_snd( jps_tice, isec, ztmp3, info )
1375         IF( ssnd(jps_tmix)%laction )   CALL cpl_prism_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1376      ENDIF
1377      !
1378      !                                                      ! ------------------------- !
1379      !                                                      !           Albedo          !
1380      !                                                      ! ------------------------- !
1381      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
1382         ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1383         CALL cpl_prism_snd( jps_albice, isec, ztmp3, info )
1384      ENDIF
1385      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
1386         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
1387         DO jl=1,jpl
1388            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
1389         ENDDO
1390         CALL cpl_prism_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1391      ENDIF
1392      !                                                      ! ------------------------- !
1393      !                                                      !  Ice fraction & Thickness !
1394      !                                                      ! ------------------------- !
1395      ! Send ice fraction field
1396      IF( ssnd(jps_fice)%laction ) THEN
1397         SELECT CASE( sn_snd_thick%clcat )
1398         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
1399         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
1400         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1401         END SELECT
1402         CALL cpl_prism_snd( jps_fice, isec, ztmp3, info )
1403      ENDIF
1404
1405      ! Send ice and snow thickness field
1406      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
1407         SELECT CASE( sn_snd_thick%cldes)
1408         CASE( 'none'                  )       ! nothing to do
1409         CASE( 'weighted ice and snow' )   
1410            SELECT CASE( sn_snd_thick%clcat )
1411            CASE( 'yes' )   
1412               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
1413               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
1414            CASE( 'no' )
1415               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
1416               DO jl=1,jpl
1417                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
1418                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
1419               ENDDO
1420            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1421            END SELECT
1422         CASE( 'ice and snow'         )   
1423            ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
1424            ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
1425         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
1426         END SELECT
1427         IF( ssnd(jps_hice)%laction )   CALL cpl_prism_snd( jps_hice, isec, ztmp3, info )
1428         IF( ssnd(jps_hsnw)%laction )   CALL cpl_prism_snd( jps_hsnw, isec, ztmp4, info )
1429      ENDIF
1430      !
1431#if defined key_cpl_carbon_cycle
1432      !                                                      ! ------------------------- !
1433      !                                                      !  CO2 flux from PISCES     !
1434      !                                                      ! ------------------------- !
1435      IF( ssnd(jps_co2)%laction )   CALL cpl_prism_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
1436      !
1437#endif
1438      !                                                      ! ------------------------- !
1439      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
1440         !                                                   ! ------------------------- !
1441         !   
1442         !                                                  j+1   j     -----V---F
1443         ! surface velocity always sent from T point                     !       |
1444         !                                                        j      |   T   U
1445         !                                                               |       |
1446         !                                                   j    j-1   -I-------|
1447         !                                               (for I)         |       |
1448         !                                                              i-1  i   i
1449         !                                                               i      i+1 (for I)
1450         SELECT CASE( TRIM( sn_snd_crt%cldes ) )
1451         CASE( 'oce only'             )      ! C-grid ==> T
1452            DO jj = 2, jpjm1
1453               DO ji = fs_2, fs_jpim1   ! vector opt.
1454                  zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
1455                  zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
1456               END DO
1457            END DO
1458         CASE( 'weighted oce and ice' )   
1459            SELECT CASE ( cp_ice_msh )
1460            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1461               DO jj = 2, jpjm1
1462                  DO ji = fs_2, fs_jpim1   ! vector opt.
1463                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1464                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
1465                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1466                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1467                  END DO
1468               END DO
1469            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1470               DO jj = 2, jpjm1
1471                  DO ji = 2, jpim1   ! NO vector opt.
1472                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1473                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1474                     zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1475                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1476                     zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1477                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1478                  END DO
1479               END DO
1480            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1481               DO jj = 2, jpjm1
1482                  DO ji = 2, jpim1   ! NO vector opt.
1483                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1484                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1485                     zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1486                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1487                     zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1488                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1489                  END DO
1490               END DO
1491            END SELECT
1492            CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
1493         CASE( 'mixed oce-ice'        )
1494            SELECT CASE ( cp_ice_msh )
1495            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1496               DO jj = 2, jpjm1
1497                  DO ji = fs_2, fs_jpim1   ! vector opt.
1498                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
1499                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1500                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
1501                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1502                  END DO
1503               END DO
1504            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1505               DO jj = 2, jpjm1
1506                  DO ji = 2, jpim1   ! NO vector opt.
1507                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1508                        &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1509                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1510                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1511                        &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1512                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1513                  END DO
1514               END DO
1515            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1516               DO jj = 2, jpjm1
1517                  DO ji = 2, jpim1   ! NO vector opt.
1518                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1519                        &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1520                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1521                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1522                        &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1523                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1524                  END DO
1525               END DO
1526            END SELECT
1527         END SELECT
1528         CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
1529         !
1530         !
1531         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
1532            !                                                                     ! Ocean component
1533            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
1534            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
1535            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
1536            zoty1(:,:) = ztmp2(:,:)
1537            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
1538               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
1539               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
1540               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
1541               zity1(:,:) = ztmp2(:,:)
1542            ENDIF
1543         ENDIF
1544         !
1545         ! spherical coordinates to cartesian -> 2 components to 3 components
1546         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
1547            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
1548            ztmp2(:,:) = zoty1(:,:)
1549            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
1550            !
1551            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
1552               ztmp1(:,:) = zitx1(:,:)
1553               ztmp1(:,:) = zity1(:,:)
1554               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
1555            ENDIF
1556         ENDIF
1557         !
1558         IF( ssnd(jps_ocx1)%laction )   CALL cpl_prism_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
1559         IF( ssnd(jps_ocy1)%laction )   CALL cpl_prism_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
1560         IF( ssnd(jps_ocz1)%laction )   CALL cpl_prism_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
1561         !
1562         IF( ssnd(jps_ivx1)%laction )   CALL cpl_prism_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
1563         IF( ssnd(jps_ivy1)%laction )   CALL cpl_prism_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
1564         IF( ssnd(jps_ivz1)%laction )   CALL cpl_prism_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
1565         !
1566      ENDIF
1567      !
1568      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1569      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1570      !
1571      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
1572      !
1573   END SUBROUTINE sbc_cpl_snd
1574   
1575#else
1576   !!----------------------------------------------------------------------
1577   !!   Dummy module                                            NO coupling
1578   !!----------------------------------------------------------------------
1579   USE par_kind        ! kind definition
1580CONTAINS
1581   SUBROUTINE sbc_cpl_snd( kt )
1582      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt
1583   END SUBROUTINE sbc_cpl_snd
1584   !
1585   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
1586      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt, k_fsbc, k_ice
1587   END SUBROUTINE sbc_cpl_rcv
1588   !
1589   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1590      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1591      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1592      p_taui(:,:) = 0.   ;   p_tauj(:,:) = 0. ! stupid definition to avoid warning message when compiling...
1593      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?'
1594   END SUBROUTINE sbc_cpl_ice_tau
1595   !
1596   SUBROUTINE sbc_cpl_ice_flx( p_frld , palbi   , psst    , pist  )
1597      REAL(wp), INTENT(in   ), DIMENSION(:,:  ) ::   p_frld     ! lead fraction                [0 to 1]
1598      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1599      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature      [Celcius]
1600      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature      [Kelvin]
1601      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', p_frld(1,1), palbi(1,1,1), psst(1,1), pist(1,1,1) 
1602   END SUBROUTINE sbc_cpl_ice_flx
1603   
1604#endif
1605
1606   !!======================================================================
1607END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.