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

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 4066

Last change on this file since 4066 was 4066, checked in by smasson, 11 years ago

trunk: bugfix for mixed oce-ice coupling, see #1110 and #1137

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