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

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

source: branches/UKMO/dev_r5518_ww3_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 6594

Last change on this file since 6594 was 6594, checked in by jcastill, 8 years ago

Fix the water level coupling by adding the bathymetry to the coupling field

File size: 135.3 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   !!----------------------------------------------------------------------
12   !!   namsbc_cpl      : coupled formulation namlist
13   !!   sbc_cpl_init    : initialisation of the coupled exchanges
14   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
15   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
16   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
17   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
18   !!   sbc_cpl_snd     : send     fields to the atmosphere
19   !!----------------------------------------------------------------------
20   USE dom_oce         ! ocean space and time domain
21   USE sbc_oce         ! Surface boundary condition: ocean fields
22   USE sbc_ice         ! Surface boundary condition: ice fields
23   USE sbcapr
24   USE sbcdcy          ! surface boundary condition: diurnal cycle
25   USE phycst          ! physical constants
26#if defined key_lim3
27   USE ice             ! ice variables
28#endif
29#if defined key_lim2
30   USE par_ice_2       ! ice parameters
31   USE ice_2           ! ice variables
32#endif
33   USE cpl_oasis3      ! OASIS3 coupling
34   USE geo2ocean       !
35   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
36   USE albedo          !
37   USE in_out_manager  ! I/O manager
38   USE iom             ! NetCDF library
39   USE lib_mpp         ! distribued memory computing library
40   USE wrk_nemo        ! work arrays
41   USE timing          ! Timing
42   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
43   USE eosbn2
44   USE sbcrnf   , ONLY : l_rnfcpl
45#if defined key_cpl_carbon_cycle
46   USE p4zflx, ONLY : oce_co2
47#endif
48#if defined key_cice
49   USE ice_domain_size, only: ncat
50#endif
51#if defined key_lim3
52   USE limthd_dh       ! for CALL lim_thd_snwblow
53#endif
54
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC   sbc_cpl_init       ! routine called by sbcmod.F90
59   PUBLIC   sbc_cpl_rcv        ! routine called by sbc_ice_lim(_2).F90
60   PUBLIC   sbc_cpl_snd        ! routine called by step.F90
61   PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
62   PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).F90
63   PUBLIC   sbc_cpl_alloc      ! routine called in sbcice_cice.F90
64
65   INTEGER, PARAMETER ::   jpr_otx1   =  1            ! 3 atmosphere-ocean stress components on grid 1
66   INTEGER, PARAMETER ::   jpr_oty1   =  2            !
67   INTEGER, PARAMETER ::   jpr_otz1   =  3            !
68   INTEGER, PARAMETER ::   jpr_otx2   =  4            ! 3 atmosphere-ocean stress components on grid 2
69   INTEGER, PARAMETER ::   jpr_oty2   =  5            !
70   INTEGER, PARAMETER ::   jpr_otz2   =  6            !
71   INTEGER, PARAMETER ::   jpr_itx1   =  7            ! 3 atmosphere-ice   stress components on grid 1
72   INTEGER, PARAMETER ::   jpr_ity1   =  8            !
73   INTEGER, PARAMETER ::   jpr_itz1   =  9            !
74   INTEGER, PARAMETER ::   jpr_itx2   = 10            ! 3 atmosphere-ice   stress components on grid 2
75   INTEGER, PARAMETER ::   jpr_ity2   = 11            !
76   INTEGER, PARAMETER ::   jpr_itz2   = 12            !
77   INTEGER, PARAMETER ::   jpr_qsroce = 13            ! Qsr above the ocean
78   INTEGER, PARAMETER ::   jpr_qsrice = 14            ! Qsr above the ice
79   INTEGER, PARAMETER ::   jpr_qsrmix = 15 
80   INTEGER, PARAMETER ::   jpr_qnsoce = 16            ! Qns above the ocean
81   INTEGER, PARAMETER ::   jpr_qnsice = 17            ! Qns above the ice
82   INTEGER, PARAMETER ::   jpr_qnsmix = 18
83   INTEGER, PARAMETER ::   jpr_rain   = 19            ! total liquid precipitation (rain)
84   INTEGER, PARAMETER ::   jpr_snow   = 20            ! solid precipitation over the ocean (snow)
85   INTEGER, PARAMETER ::   jpr_tevp   = 21            ! total evaporation
86   INTEGER, PARAMETER ::   jpr_ievp   = 22            ! solid evaporation (sublimation)
87   INTEGER, PARAMETER ::   jpr_sbpr   = 23            ! sublimation - liquid precipitation - solid precipitation
88   INTEGER, PARAMETER ::   jpr_semp   = 24            ! solid freshwater budget (sublimation - snow)
89   INTEGER, PARAMETER ::   jpr_oemp   = 25            ! ocean freshwater budget (evap - precip)
90   INTEGER, PARAMETER ::   jpr_w10m   = 26            ! 10m wind
91   INTEGER, PARAMETER ::   jpr_dqnsdt = 27            ! d(Q non solar)/d(temperature)
92   INTEGER, PARAMETER ::   jpr_rnf    = 28            ! runoffs
93   INTEGER, PARAMETER ::   jpr_cal    = 29            ! calving
94   INTEGER, PARAMETER ::   jpr_taum   = 30            ! wind stress module
95   INTEGER, PARAMETER ::   jpr_co2    = 31
96   INTEGER, PARAMETER ::   jpr_topm   = 32            ! topmeltn
97   INTEGER, PARAMETER ::   jpr_botm   = 33            ! botmeltn
98   INTEGER, PARAMETER ::   jpr_sflx   = 34            ! salt flux
99   INTEGER, PARAMETER ::   jpr_toce   = 35            ! ocean temperature
100   INTEGER, PARAMETER ::   jpr_soce   = 36            ! ocean salinity
101   INTEGER, PARAMETER ::   jpr_ocx1   = 37            ! ocean current on grid 1
102   INTEGER, PARAMETER ::   jpr_ocy1   = 38            !
103   INTEGER, PARAMETER ::   jpr_ssh    = 39            ! sea surface height
104   INTEGER, PARAMETER ::   jpr_fice   = 40            ! ice fraction         
105   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness
106   INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level
107   INTEGER, PARAMETER ::   jpr_mslp   = 43            ! mean sea level pressure
108   INTEGER, PARAMETER ::   jpr_hsig   = 44            ! Hsig
109   INTEGER, PARAMETER ::   jpr_phioc  = 45            ! Wave=>ocean energy flux
110   INTEGER, PARAMETER ::   jprcv      = 45            ! total number of fields received
111
112   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
113   INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
114   INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
115   INTEGER, PARAMETER ::   jps_tmix   =  4            ! mixed temperature (ocean+ice)
116   INTEGER, PARAMETER ::   jps_albice =  5            ! ice   albedo
117   INTEGER, PARAMETER ::   jps_albmix =  6            ! mixed albedo
118   INTEGER, PARAMETER ::   jps_hice   =  7            ! ice  thickness
119   INTEGER, PARAMETER ::   jps_hsnw   =  8            ! snow thickness
120   INTEGER, PARAMETER ::   jps_ocx1   =  9            ! ocean current on grid 1
121   INTEGER, PARAMETER ::   jps_ocy1   = 10            !
122   INTEGER, PARAMETER ::   jps_ocz1   = 11            !
123   INTEGER, PARAMETER ::   jps_ivx1   = 12            ! ice   current on grid 1
124   INTEGER, PARAMETER ::   jps_ivy1   = 13            !
125   INTEGER, PARAMETER ::   jps_ivz1   = 14            !
126   INTEGER, PARAMETER ::   jps_co2    = 15
127   INTEGER, PARAMETER ::   jps_soce   = 16            ! ocean salinity
128   INTEGER, PARAMETER ::   jps_ssh    = 17            ! sea surface height
129   INTEGER, PARAMETER ::   jps_qsroce = 18            ! Qsr above the ocean
130   INTEGER, PARAMETER ::   jps_qnsoce = 19            ! Qns above the ocean
131   INTEGER, PARAMETER ::   jps_oemp   = 20            ! ocean freshwater budget (evap - precip)
132   INTEGER, PARAMETER ::   jps_sflx   = 21            ! salt flux
133   INTEGER, PARAMETER ::   jps_otx1   = 22            ! 2 atmosphere-ocean stress components on grid 1
134   INTEGER, PARAMETER ::   jps_oty1   = 23            !
135   INTEGER, PARAMETER ::   jps_rnf    = 24            ! runoffs
136   INTEGER, PARAMETER ::   jps_taum   = 25            ! wind stress module
137   INTEGER, PARAMETER ::   jps_fice2  = 26            ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling)
138   INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl)
139   INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level
140   INTEGER, PARAMETER ::   jps_ficet  = 29            ! total ice fraction
141   INTEGER, PARAMETER ::   jps_ocxw   = 30            ! currents on grid 1
142   INTEGER, PARAMETER ::   jps_ocyw   = 31            ! currents on grid 1
143   INTEGER, PARAMETER ::   jps_wlev   = 32            ! water level
144   INTEGER, PARAMETER ::   jpsnd      = 32            ! total number of fields sended
145
146   !                                                         !!** namelist namsbc_cpl **
147   TYPE ::   FLD_C
148      CHARACTER(len = 32) ::   cldes                  ! desciption of the coupling strategy
149      CHARACTER(len = 32) ::   clcat                  ! multiple ice categories strategy
150      CHARACTER(len = 32) ::   clvref                 ! reference of vector ('spherical' or 'cartesian')
151      CHARACTER(len = 32) ::   clvor                  ! orientation of vector fields ('eastward-northward' or 'local grid')
152      CHARACTER(len = 32) ::   clvgrd                 ! grids on which is located the vector fields
153   END TYPE FLD_C
154   ! Send to the atmosphere                           !
155   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                       
156   ! Received from the atmosphere                     !
157   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
158   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp                         
159   ! Send to waves
160   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev
161   ! Received from waves
162   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc
163   ! Other namelist parameters                        !
164   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
165   LOGICAL     ::   ln_usecplmask          !  use a coupling mask file to merge data received from several models
166                                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
167   TYPE ::   DYNARR     
168      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
169   END TYPE DYNARR
170
171   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                      ! all fields recieved from the atmosphere
172
173   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
174
175   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]
176   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0)
177
178   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
179
180   !! Substitution
181#  include "domzgr_substitute.h90"
182#  include "vectopt_loop_substitute.h90"
183   !!----------------------------------------------------------------------
184   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
185   !! $Id$
186   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
187   !!----------------------------------------------------------------------
188
189CONTAINS
190 
191   INTEGER FUNCTION sbc_cpl_alloc()
192      !!----------------------------------------------------------------------
193      !!             ***  FUNCTION sbc_cpl_alloc  ***
194      !!----------------------------------------------------------------------
195      INTEGER :: ierr(4)
196      !!----------------------------------------------------------------------
197      ierr(:) = 0
198      !
199      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
200     
201#if ! defined key_lim3 && ! defined key_lim2 && ! defined key_cice
202      ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
203#endif
204      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
205      !
206      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) )
207
208      sbc_cpl_alloc = MAXVAL( ierr )
209      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
210      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
211      !
212   END FUNCTION sbc_cpl_alloc
213
214
215   SUBROUTINE sbc_cpl_init( k_ice )     
216      !!----------------------------------------------------------------------
217      !!             ***  ROUTINE sbc_cpl_init  ***
218      !!
219      !! ** Purpose :   Initialisation of send and received information from
220      !!                the atmospheric component
221      !!
222      !! ** Method  : * Read namsbc_cpl namelist
223      !!              * define the receive interface
224      !!              * define the send    interface
225      !!              * initialise the OASIS coupler
226      !!----------------------------------------------------------------------
227      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
228      !!
229      INTEGER ::   jn   ! dummy loop index
230      INTEGER ::   ios  ! Local integer output status for namelist read
231      INTEGER ::   inum 
232      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
233      !!
234      NAMELIST/namsbc_cpl/  sn_snd_temp , sn_snd_alb  , sn_snd_thick , sn_snd_crt   , sn_snd_co2,      &
235         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      &
236         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   &
237         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,   &
238         &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp
239      !!---------------------------------------------------------------------
240      !
241      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
242      !
243      CALL wrk_alloc( jpi,jpj, zacs, zaos )
244
245      ! ================================ !
246      !      Namelist informations       !
247      ! ================================ !
248
249      REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
250      READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
251901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
252
253      REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
254      READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
255902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
256      IF(lwm) WRITE ( numond, namsbc_cpl )
257
258      IF(lwp) THEN                        ! control print
259         WRITE(numout,*)
260         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
261         WRITE(numout,*)'~~~~~~~~~~~~'
262      ENDIF
263      IF( lwp .AND. ln_cpl ) THEN                        ! control print
264         WRITE(numout,*)'  received fields (mutiple ice categogies)'
265         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
266         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
267         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
268         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
269         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
270         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
271         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
272         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
273         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
274         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
275         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
276         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
277         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
278         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
279         WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'
280         WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'
281         WRITE(numout,*)'  sent fields (multiple ice categories)'
282         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
283         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
284         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
285         WRITE(numout,*)'      total ice fraction              = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%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         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')'
292         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')'
293         WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes   ), ' (', TRIM(sn_snd_crtw%clcat   ), ')'
294         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref
295         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor
296         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd
297         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
298         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
299      ENDIF
300
301      !                                   ! allocate sbccpl arrays
302      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
303     
304      ! ================================ !
305      !   Define the receive interface   !
306      ! ================================ !
307      nrcvinfo(:) = OASIS_idle   ! needed by nrcvinfo(jpr_otx1) if we do not receive ocean stress
308
309      ! for each field: define the OASIS name                              (srcv(:)%clname)
310      !                 define receive or not from the namelist parameters (srcv(:)%laction)
311      !                 define the north fold type of lbc                  (srcv(:)%nsgn)
312
313      ! default definitions of srcv
314      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
315
316      !                                                      ! ------------------------- !
317      !                                                      ! ice and ocean wind stress !   
318      !                                                      ! ------------------------- !
319      !                                                           ! Name
320      srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
321      srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
322      srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
323      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
324      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
325      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
326      !
327      srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
328      srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
329      srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
330      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
331      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
332      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
333      !
334      ! Vectors: change of sign at north fold ONLY if on the local grid
335      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
336     
337      !                                                           ! Set grid and action
338      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
339      CASE( 'T' ) 
340         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
341         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
342         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
343      CASE( 'U,V' ) 
344         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
345         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
346         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
347         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
348         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
349      CASE( 'U,V,T' )
350         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
351         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
352         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
353         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
354         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
355      CASE( 'U,V,I' )
356         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
357         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
358         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
359         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
360         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
361      CASE( 'U,V,F' )
362         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
363         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
364         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
365         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
366         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
367      CASE( 'T,I' ) 
368         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
369         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
370         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
371         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
372      CASE( 'T,F' ) 
373         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
374         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
375         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
376         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
377      CASE( 'T,U,V' )
378         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
379         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
380         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
381         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
382         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
383      CASE default   
384         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
385      END SELECT
386      !
387      IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
388         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE. 
389      !
390      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
391            srcv(jpr_otx2:jpr_otz2)%laction = .FALSE. 
392            srcv(jpr_itx2:jpr_itz2)%laction = .FALSE. 
393            srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid   ! not needed but cleaner...
394            srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid   ! not needed but cleaner...
395      ENDIF
396      !
397      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
398         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
399         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
400         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
401      ENDIF
402       
403      !                                                      ! ------------------------- !
404      !                                                      !    freshwater budget      !   E-P
405      !                                                      ! ------------------------- !
406      ! we suppose that atmosphere modele do not make the difference between precipiration (liquide or solid)
407      ! over ice of free ocean within the same atmospheric cell.cd
408      srcv(jpr_rain)%clname = 'OTotRain'      ! Rain = liquid precipitation
409      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
410      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
411      srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation
412      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation
413      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
414      srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
415      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
416      CASE( 'none'          )       ! nothing to do
417      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
418      CASE( 'conservative'  )
419         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
420         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
421      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
422      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
423      END SELECT
424
425      !                                                      ! ------------------------- !
426      !                                                      !     Runoffs & Calving     !   
427      !                                                      ! ------------------------- !
428      srcv(jpr_rnf   )%clname = 'O_Runoff'
429      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
430         srcv(jpr_rnf)%laction = .TRUE.
431         l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
432         ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
433         IF(lwp) WRITE(numout,*)
434         IF(lwp) WRITE(numout,*) '   runoffs received from oasis -> force ln_rnf = ', ln_rnf
435      ENDIF
436      !
437      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
438
439      !                                                      ! ------------------------- !
440      !                                                      !    non solar radiation    !   Qns
441      !                                                      ! ------------------------- !
442      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
443      srcv(jpr_qnsice)%clname = 'O_QnsIce'
444      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
445      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
446      CASE( 'none'          )       ! nothing to do
447      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
448      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
449      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
450      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
451      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
452      END SELECT
453      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
454         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
455      !                                                      ! ------------------------- !
456      !                                                      !    solar radiation        !   Qsr
457      !                                                      ! ------------------------- !
458      srcv(jpr_qsroce)%clname = 'O_QsrOce'
459      srcv(jpr_qsrice)%clname = 'O_QsrIce'
460      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
461      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
462      CASE( 'none'          )       ! nothing to do
463      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
464      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
465      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
466      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
467      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
468      END SELECT
469      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
470         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
471      !                                                      ! ------------------------- !
472      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
473      !                                                      ! ------------------------- !
474      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
475      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
476      !
477      ! non solar sensitivity mandatory for LIM ice model
478      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4 .AND. nn_components /= jp_iam_sas ) &
479         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
480      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
481      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
482         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
483      !                                                      ! ------------------------- !
484      !                                                      !      10m wind module      !   
485      !                                                      ! ------------------------- !
486      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
487      !
488      !                                                      ! ------------------------- !
489      !                                                      !   wind stress module      !   
490      !                                                      ! ------------------------- !
491      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
492      lhftau = srcv(jpr_taum)%laction
493
494      !                                                      ! ------------------------- !
495      !                                                      !      Atmospheric CO2      !
496      !                                                      ! ------------------------- !
497      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
498
499      !                                                      ! ------------------------- !
500      !                                                      ! Mean Sea Level Pressure   !
501      !                                                      ! ------------------------- !
502      srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE.
503
504      !                                                      ! ------------------------- !
505      !                                                      !   topmelt and botmelt     !   
506      !                                                      ! ------------------------- !
507      srcv(jpr_topm )%clname = 'OTopMlt'
508      srcv(jpr_botm )%clname = 'OBotMlt'
509      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
510         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
511            srcv(jpr_topm:jpr_botm)%nct = jpl
512         ELSE
513            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
514         ENDIF
515         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
516      ENDIF
517      !                                                      ! ------------------------- !                                                      !
518      !                                                      !      Wave breaking        !   
519      !                                                      ! ------------------------- !
520      srcv(jpr_hsig)%clname  = 'O_Hsigwa'  ;   IF( TRIM(sn_rcv_hsig%cldes  ) == 'coupled' )   srcv(jpr_hsig)%laction = .TRUE.
521      srcv(jpr_phioc)%clname = 'O_PhiOce'  ;   IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' )   srcv(jpr_phioc)%laction = .TRUE.
522      !
523      !                                                      ! ------------------------------- !
524      !                                                      !   OPA-SAS coupling - rcv by opa !   
525      !                                                      ! ------------------------------- !
526      srcv(jpr_sflx)%clname = 'O_SFLX'
527      srcv(jpr_fice)%clname = 'RIceFrc'
528      !
529      IF( nn_components == jp_iam_opa ) THEN    ! OPA coupled to SAS via OASIS: force received field by OPA (sent by SAS)
530         srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
531         srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
532         srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
533         srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
534         srcv(jpr_otx1)%clgrid = 'U'        ! oce components given at U-point
535         srcv(jpr_oty1)%clgrid = 'V'        !           and           V-point
536         ! Vectors: change of sign at north fold ONLY if on the local grid
537         srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
538         sn_rcv_tau%clvgrd = 'U,V'
539         sn_rcv_tau%clvor = 'local grid'
540         sn_rcv_tau%clvref = 'spherical'
541         sn_rcv_emp%cldes = 'oce only'
542         !
543         IF(lwp) THEN                        ! control print
544            WRITE(numout,*)
545            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
546            WRITE(numout,*)'               OPA component  '
547            WRITE(numout,*)
548            WRITE(numout,*)'  received fields from SAS component '
549            WRITE(numout,*)'                  ice cover '
550            WRITE(numout,*)'                  oce only EMP  '
551            WRITE(numout,*)'                  salt flux  '
552            WRITE(numout,*)'                  mixed oce-ice solar flux  '
553            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
554            WRITE(numout,*)'                  wind stress U,V on local grid and sperical coordinates '
555            WRITE(numout,*)'                  wind stress module'
556            WRITE(numout,*)
557         ENDIF
558      ENDIF
559      !                                                      ! -------------------------------- !
560      !                                                      !   OPA-SAS coupling - rcv by sas  !   
561      !                                                      ! -------------------------------- !
562      srcv(jpr_toce  )%clname = 'I_SSTSST'
563      srcv(jpr_soce  )%clname = 'I_SSSal'
564      srcv(jpr_ocx1  )%clname = 'I_OCurx1'
565      srcv(jpr_ocy1  )%clname = 'I_OCury1'
566      srcv(jpr_ssh   )%clname = 'I_SSHght'
567      srcv(jpr_e3t1st)%clname = 'I_E3T1st'   
568      srcv(jpr_fraqsr)%clname = 'I_FraQsr'   
569      !
570      IF( nn_components == jp_iam_sas ) THEN
571         IF( .NOT. ln_cpl ) srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
572         IF( .NOT. ln_cpl ) srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
573         IF( .NOT. ln_cpl ) srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
574         srcv( (/jpr_toce, jpr_soce, jpr_ssh, jpr_fraqsr, jpr_ocx1, jpr_ocy1/) )%laction = .TRUE.
575         srcv( jpr_e3t1st )%laction = lk_vvl
576         srcv(jpr_ocx1)%clgrid = 'U'        ! oce components given at U-point
577         srcv(jpr_ocy1)%clgrid = 'V'        !           and           V-point
578         ! Vectors: change of sign at north fold ONLY if on the local grid
579         srcv(jpr_ocx1:jpr_ocy1)%nsgn = -1.
580         ! Change first letter to couple with atmosphere if already coupled OPA
581         ! this is nedeed as each variable name used in the namcouple must be unique:
582         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
583         DO jn = 1, jprcv
584            IF ( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
585         END DO
586         !
587         IF(lwp) THEN                        ! control print
588            WRITE(numout,*)
589            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
590            WRITE(numout,*)'               SAS component  '
591            WRITE(numout,*)
592            IF( .NOT. ln_cpl ) THEN
593               WRITE(numout,*)'  received fields from OPA component '
594            ELSE
595               WRITE(numout,*)'  Additional received fields from OPA component : '
596            ENDIF
597            WRITE(numout,*)'               sea surface temperature (Celcius) '
598            WRITE(numout,*)'               sea surface salinity ' 
599            WRITE(numout,*)'               surface currents ' 
600            WRITE(numout,*)'               sea surface height ' 
601            WRITE(numout,*)'               thickness of first ocean T level '       
602            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
603            WRITE(numout,*)
604         ENDIF
605      ENDIF
606     
607      ! =================================================== !
608      ! Allocate all parts of frcv used for received fields !
609      ! =================================================== !
610      DO jn = 1, jprcv
611         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
612      END DO
613      ! Allocate taum part of frcv which is used even when not received as coupling field
614      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
615      ! Allocate w10m part of frcv which is used even when not received as coupling field
616      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
617      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
618      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
619      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
620      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
621      IF( k_ice /= 0 ) THEN
622         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
623         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
624      END IF
625
626      ! ================================ !
627      !     Define the send interface    !
628      ! ================================ !
629      ! for each field: define the OASIS name                           (ssnd(:)%clname)
630      !                 define send or not from the namelist parameters (ssnd(:)%laction)
631      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
632     
633      ! default definitions of nsnd
634      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
635         
636      !                                                      ! ------------------------- !
637      !                                                      !    Surface temperature    !
638      !                                                      ! ------------------------- !
639      ssnd(jps_toce)%clname = 'O_SSTSST'
640      ssnd(jps_tice)%clname = 'O_TepIce'
641      ssnd(jps_tmix)%clname = 'O_TepMix'
642      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
643      CASE( 'none'                                 )       ! nothing to do
644      CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE.
645      CASE( 'oce and ice' , 'weighted oce and ice' )
646         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
647         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
648      CASE( 'mixed oce-ice'                        )   ;   ssnd( jps_tmix )%laction = .TRUE.
649      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
650      END SELECT
651           
652      !                                                      ! ------------------------- !
653      !                                                      !          Albedo           !
654      !                                                      ! ------------------------- !
655      ssnd(jps_albice)%clname = 'O_AlbIce' 
656      ssnd(jps_albmix)%clname = 'O_AlbMix'
657      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
658      CASE( 'none'                 )     ! nothing to do
659      CASE( 'ice' , 'weighted ice' )   ; ssnd(jps_albice)%laction = .TRUE.
660      CASE( 'mixed oce-ice'        )   ; ssnd(jps_albmix)%laction = .TRUE.
661      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
662      END SELECT
663      !
664      ! Need to calculate oceanic albedo if
665      !     1. sending mixed oce-ice albedo or
666      !     2. receiving mixed oce-ice solar radiation
667      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
668         CALL albedo_oce( zaos, zacs )
669         ! Due to lack of information on nebulosity : mean clear/overcast sky
670         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
671      ENDIF
672
673      !                                                      ! ------------------------- !
674      !                                                      !  Ice fraction & Thickness !
675      !                                                      ! ------------------------- !
676      ssnd(jps_fice)%clname = 'OIceFrc'
677      ssnd(jps_ficet)%clname = 'OIceFrcT'
678      ssnd(jps_hice)%clname = 'OIceTck'
679      ssnd(jps_hsnw)%clname = 'OSnwTck'
680      IF( k_ice /= 0 ) THEN
681         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
682! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
683         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
684      ENDIF
685     
686      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled')  ssnd(jps_ficet)%laction = .TRUE.
687
688      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
689      CASE( 'none'         )       ! nothing to do
690      CASE( 'ice and snow' ) 
691         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
692         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
693            ssnd(jps_hice:jps_hsnw)%nct = jpl
694         ENDIF
695      CASE ( 'weighted ice and snow' ) 
696         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
697         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
698      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
699      END SELECT
700
701      !                                                      ! ------------------------- !
702      !                                                      !      Surface current      !
703      !                                                      ! ------------------------- !
704      !        ocean currents              !            ice velocities
705      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
706      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
707      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
708      ssnd(jps_ocxw)%clname = 'O_OCurxw'
709      ssnd(jps_ocyw)%clname = 'O_OCuryw'
710      !
711      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
712
713      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
714         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
715      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
716         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
717         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
718      ENDIF
719      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
720      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
721      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
722      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
723      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
724      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
725      CASE( 'weighted oce and ice' )   !   nothing to do
726      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
727      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
728      END SELECT
729
730      ssnd(jps_ocxw:jps_ocyw)%nsgn = -1.   ! vectors: change of the sign at the north fold
731
732      IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN
733         ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V'
734      ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN
735         CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' )
736      ENDIF
737      IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1.
738      SELECT CASE( TRIM( sn_snd_crtw%cldes ) )
739      CASE( 'none'                 )   ;   ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE.
740      CASE( 'oce only'             )   ;   ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE.
741      CASE( 'weighted oce and ice' )   !   nothing to do
742      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
743      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' )
744      END SELECT
745
746      !                                                      ! ------------------------- !
747      !                                                      !          CO2 flux         !
748      !                                                      ! ------------------------- !
749      ssnd(jps_co2)%clname  = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
750
751      !                                                      ! ------------------------- !
752      !                                                      !     Sea surface height    !
753      !                                                      ! ------------------------- !
754      ssnd(jps_wlev)%clname = 'O_Wlevel' ;  IF( TRIM(sn_snd_wlev%cldes) == 'coupled' )   ssnd(jps_wlev)%laction = .TRUE.
755
756      !                                                      ! ------------------------------- !
757      !                                                      !   OPA-SAS coupling - snd by opa !   
758      !                                                      ! ------------------------------- !
759      ssnd(jps_ssh   )%clname = 'O_SSHght'
760      ssnd(jps_soce  )%clname = 'O_SSSal' 
761      ssnd(jps_e3t1st)%clname = 'O_E3T1st'   
762      ssnd(jps_fraqsr)%clname = 'O_FraQsr'
763      !
764      IF( nn_components == jp_iam_opa ) THEN
765         ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
766         ssnd( (/jps_toce, jps_soce, jps_ssh, jps_fraqsr, jps_ocx1, jps_ocy1/) )%laction = .TRUE.
767         ssnd( jps_e3t1st )%laction = lk_vvl
768         ! vector definition: not used but cleaner...
769         ssnd(jps_ocx1)%clgrid  = 'U'        ! oce components given at U-point
770         ssnd(jps_ocy1)%clgrid  = 'V'        !           and           V-point
771         sn_snd_crt%clvgrd = 'U,V'
772         sn_snd_crt%clvor = 'local grid'
773         sn_snd_crt%clvref = 'spherical'
774         !
775         IF(lwp) THEN                        ! control print
776            WRITE(numout,*)
777            WRITE(numout,*)'  sent fields to SAS component '
778            WRITE(numout,*)'               sea surface temperature (T before, Celcius) '
779            WRITE(numout,*)'               sea surface salinity ' 
780            WRITE(numout,*)'               surface currents U,V on local grid and spherical coordinates' 
781            WRITE(numout,*)'               sea surface height ' 
782            WRITE(numout,*)'               thickness of first ocean T level '       
783            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
784            WRITE(numout,*)
785         ENDIF
786      ENDIF
787      !                                                      ! ------------------------------- !
788      !                                                      !   OPA-SAS coupling - snd by sas !   
789      !                                                      ! ------------------------------- !
790      ssnd(jps_sflx  )%clname = 'I_SFLX'     
791      ssnd(jps_fice2 )%clname = 'IIceFrc'
792      ssnd(jps_qsroce)%clname = 'I_QsrOce'   
793      ssnd(jps_qnsoce)%clname = 'I_QnsOce'   
794      ssnd(jps_oemp  )%clname = 'IOEvaMPr' 
795      ssnd(jps_otx1  )%clname = 'I_OTaux1'   
796      ssnd(jps_oty1  )%clname = 'I_OTauy1'   
797      ssnd(jps_rnf   )%clname = 'I_Runoff'   
798      ssnd(jps_taum  )%clname = 'I_TauMod'   
799      !
800      IF( nn_components == jp_iam_sas ) THEN
801         IF( .NOT. ln_cpl ) ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
802         ssnd( (/jps_qsroce, jps_qnsoce, jps_oemp, jps_fice2, jps_sflx, jps_otx1, jps_oty1, jps_taum/) )%laction = .TRUE.
803         !
804         ! Change first letter to couple with atmosphere if already coupled with sea_ice
805         ! this is nedeed as each variable name used in the namcouple must be unique:
806         ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere
807         DO jn = 1, jpsnd
808            IF ( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname))
809         END DO
810         !
811         IF(lwp) THEN                        ! control print
812            WRITE(numout,*)
813            IF( .NOT. ln_cpl ) THEN
814               WRITE(numout,*)'  sent fields to OPA component '
815            ELSE
816               WRITE(numout,*)'  Additional sent fields to OPA component : '
817            ENDIF
818            WRITE(numout,*)'                  ice cover '
819            WRITE(numout,*)'                  oce only EMP  '
820            WRITE(numout,*)'                  salt flux  '
821            WRITE(numout,*)'                  mixed oce-ice solar flux  '
822            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
823            WRITE(numout,*)'                  wind stress U,V components'
824            WRITE(numout,*)'                  wind stress module'
825         ENDIF
826      ENDIF
827
828      !
829      ! ================================ !
830      !   initialisation of the coupler  !
831      ! ================================ !
832
833      CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
834     
835      IF (ln_usecplmask) THEN
836         xcplmask(:,:,:) = 0.
837         CALL iom_open( 'cplmask', inum )
838         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:nlci,1:nlcj,1:nn_cplmodel),   &
839            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nn_cplmodel /) )
840         CALL iom_close( inum )
841      ELSE
842         xcplmask(:,:,:) = 1.
843      ENDIF
844      xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
845      !
846      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
847      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
848         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
849      ncpl_qsr_freq = 86400 / ncpl_qsr_freq
850
851      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
852      !
853      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
854      !
855   END SUBROUTINE sbc_cpl_init
856
857
858   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
859      !!----------------------------------------------------------------------
860      !!             ***  ROUTINE sbc_cpl_rcv  ***
861      !!
862      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
863      !!                provide the ocean heat and freshwater fluxes.
864      !!
865      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
866      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
867      !!                to know if the field was really received or not
868      !!
869      !!              --> If ocean stress was really received:
870      !!
871      !!                  - transform the received ocean stress vector from the received
872      !!                 referential and grid into an atmosphere-ocean stress in
873      !!                 the (i,j) ocean referencial and at the ocean velocity point.
874      !!                    The received stress are :
875      !!                     - defined by 3 components (if cartesian coordinate)
876      !!                            or by 2 components (if spherical)
877      !!                     - oriented along geographical   coordinate (if eastward-northward)
878      !!                            or  along the local grid coordinate (if local grid)
879      !!                     - given at U- and V-point, resp.   if received on 2 grids
880      !!                            or at T-point               if received on 1 grid
881      !!                    Therefore and if necessary, they are successively
882      !!                  processed in order to obtain them
883      !!                     first  as  2 components on the sphere
884      !!                     second as  2 components oriented along the local grid
885      !!                     third  as  2 components on the U,V grid
886      !!
887      !!              -->
888      !!
889      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
890      !!             and total ocean freshwater fluxes 
891      !!
892      !! ** Method  :   receive all fields from the atmosphere and transform
893      !!              them into ocean surface boundary condition fields
894      !!
895      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
896      !!                        taum         wind stress module at T-point
897      !!                        wndm         wind speed  module at T-point over free ocean or leads in presence of sea-ice
898      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
899      !!                                     and the latent heat flux of solid precip. melting
900      !!                        qsr          solar ocean heat fluxes   (ocean only case)
901      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
902      !!----------------------------------------------------------------------
903      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
904      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation
905      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
906
907      !!
908      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
909      INTEGER  ::   ji, jj, jn             ! dummy loop indices
910      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
911      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
912      REAL(wp) ::   zcoef                  ! temporary scalar
913      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
914      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
915      REAL(wp) ::   zzx, zzy               ! temporary variables
916      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr
917      !!----------------------------------------------------------------------
918      !
919      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
920      !
921      CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
922      !
923      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
924      !
925      !                                                      ! ======================================================= !
926      !                                                      ! Receive all the atmos. fields (including ice information)
927      !                                                      ! ======================================================= !
928      isec = ( kt - nit000 ) * NINT( rdttra(1) )                ! date of exchanges
929      DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
930         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
931      END DO
932
933      !                                                      ! ========================= !
934      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
935         !                                                   ! ========================= !
936         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
937         ! => need to be done only when we receive the field
938         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
939            !
940            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
941               !                                                       ! (cartesian to spherical -> 3 to 2 components)
942               !
943               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
944                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
945               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
946               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
947               !
948               IF( srcv(jpr_otx2)%laction ) THEN
949                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
950                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
951                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
952                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
953               ENDIF
954               !
955            ENDIF
956            !
957            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
958               !                                                       ! (geographical to local grid -> rotate the components)
959               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
960               IF( srcv(jpr_otx2)%laction ) THEN
961                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
962               ELSE 
963                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
964               ENDIF
965               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
966               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
967            ENDIF
968            !                             
969            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
970               DO jj = 2, jpjm1                                          ! T ==> (U,V)
971                  DO ji = fs_2, fs_jpim1   ! vector opt.
972                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
973                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
974                  END DO
975               END DO
976               CALL lbc_lnk( frcv(jpr_otx1)%z3(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
977            ENDIF
978            llnewtx = .TRUE.
979         ELSE
980            llnewtx = .FALSE.
981         ENDIF
982         !                                                   ! ========================= !
983      ELSE                                                   !   No dynamical coupling   !
984         !                                                   ! ========================= !
985         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
986         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
987         llnewtx = .TRUE.
988         !
989      ENDIF
990      !                                                      ! ========================= !
991      !                                                      !    wind stress module     !   (taum)
992      !                                                      ! ========================= !
993      !
994      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
995         ! => need to be done only when otx1 was changed
996         IF( llnewtx ) THEN
997!CDIR NOVERRCHK
998            DO jj = 2, jpjm1
999!CDIR NOVERRCHK
1000               DO ji = fs_2, fs_jpim1   ! vect. opt.
1001                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
1002                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
1003                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
1004               END DO
1005            END DO
1006            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
1007            llnewtau = .TRUE.
1008         ELSE
1009            llnewtau = .FALSE.
1010         ENDIF
1011      ELSE
1012         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
1013         ! Stress module can be negative when received (interpolation problem)
1014         IF( llnewtau ) THEN
1015            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
1016         ENDIF
1017      ENDIF
1018      !
1019      !                                                      ! ========================= !
1020      !                                                      !      10 m wind speed      !   (wndm)
1021      !                                                      ! ========================= !
1022      !
1023      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
1024         ! => need to be done only when taumod was changed
1025         IF( llnewtau ) THEN
1026            zcoef = 1. / ( zrhoa * zcdrag ) 
1027!CDIR NOVERRCHK
1028            DO jj = 1, jpj
1029!CDIR NOVERRCHK
1030               DO ji = 1, jpi 
1031                  frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
1032               END DO
1033            END DO
1034         ENDIF
1035      ENDIF
1036
1037      ! u(v)tau and taum will be modified by ice model
1038      ! -> need to be reset before each call of the ice/fsbc     
1039      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
1040         !
1041         IF( ln_mixcpl ) THEN
1042            utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)
1043            vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
1044            taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
1045            wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
1046         ELSE
1047            utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
1048            vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
1049            taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
1050            wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
1051         ENDIF
1052         CALL iom_put( "taum_oce", taum )   ! output wind stress module
1053         
1054      ENDIF
1055
1056#if defined key_cpl_carbon_cycle
1057      !                                                      ! ================== !
1058      !                                                      ! atmosph. CO2 (ppm) !
1059      !                                                      ! ================== !
1060      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
1061#endif
1062      !
1063      !                                                      ! ========================= !
1064      !                                                      ! Mean Sea Level Pressure   !   (taum)
1065      !                                                      ! ========================= !
1066      !
1067      IF( srcv(jpr_mslp)%laction ) THEN                    ! UKMO SHELF effect of atmospheric pressure on SSH
1068         IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields
1069
1070         r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization
1071         ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer)
1072         apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1)                        !atmospheric pressure
1073
1074         IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible)
1075      END IF
1076
1077      !  Fields received by SAS when OASIS coupling
1078      !  (arrays no more filled at sbcssm stage)
1079      !                                                      ! ================== !
1080      !                                                      !        SSS         !
1081      !                                                      ! ================== !
1082      IF( srcv(jpr_soce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1083         sss_m(:,:) = frcv(jpr_soce)%z3(:,:,1)
1084         CALL iom_put( 'sss_m', sss_m )
1085      ENDIF
1086      !                                               
1087      !                                                      ! ================== !
1088      !                                                      !        SST         !
1089      !                                                      ! ================== !
1090      IF( srcv(jpr_toce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1091         sst_m(:,:) = frcv(jpr_toce)%z3(:,:,1)
1092         IF( srcv(jpr_soce)%laction .AND. ln_useCT ) THEN    ! make sure that sst_m is the potential temperature
1093            sst_m(:,:) = eos_pt_from_ct( sst_m(:,:), sss_m(:,:) )
1094         ENDIF
1095      ENDIF
1096      !                                                      ! ================== !
1097      !                                                      !        SSH         !
1098      !                                                      ! ================== !
1099      IF( srcv(jpr_ssh )%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1100         ssh_m(:,:) = frcv(jpr_ssh )%z3(:,:,1)
1101         CALL iom_put( 'ssh_m', ssh_m )
1102      ENDIF
1103      !                                                      ! ================== !
1104      !                                                      !  surface currents  !
1105      !                                                      ! ================== !
1106      IF( srcv(jpr_ocx1)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1107         ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
1108         ub (:,:,1) = ssu_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
1109         CALL iom_put( 'ssu_m', ssu_m )
1110      ENDIF
1111      IF( srcv(jpr_ocy1)%laction ) THEN
1112         ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
1113         vb (:,:,1) = ssv_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
1114         CALL iom_put( 'ssv_m', ssv_m )
1115      ENDIF
1116      !                                                      ! ======================== !
1117      !                                                      !  first T level thickness !
1118      !                                                      ! ======================== !
1119      IF( srcv(jpr_e3t1st )%laction ) THEN                   ! received by sas in case of opa <-> sas coupling
1120         e3t_m(:,:) = frcv(jpr_e3t1st )%z3(:,:,1)
1121         CALL iom_put( 'e3t_m', e3t_m(:,:) )
1122      ENDIF
1123      !                                                      ! ================================ !
1124      !                                                      !  fraction of solar net radiation !
1125      !                                                      ! ================================ !
1126      IF( srcv(jpr_fraqsr)%laction ) THEN                    ! received by sas in case of opa <-> sas coupling
1127         frq_m(:,:) = frcv(jpr_fraqsr)%z3(:,:,1)
1128         CALL iom_put( 'frq_m', frq_m )
1129      ENDIF
1130     
1131      !                                                      ! ========================= !
1132      IF( k_ice <= 1 .AND. MOD( kt-1, k_fsbc ) == 0 ) THEN   !  heat & freshwater fluxes ! (Ocean only case)
1133         !                                                   ! ========================= !
1134         !
1135         !                                                       ! total freshwater fluxes over the ocean (emp)
1136         IF( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction ) THEN
1137            SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
1138            CASE( 'conservative' )
1139               zemp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
1140            CASE( 'oce only', 'oce and ice' )
1141               zemp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
1142            CASE default
1143               CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
1144            END SELECT
1145         ELSE
1146            zemp(:,:) = 0._wp
1147         ENDIF
1148         !
1149         !                                                        ! runoffs and calving (added in emp)
1150         IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1151         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
1152         
1153         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
1154         ELSE                   ;   emp(:,:) =                              zemp(:,:)
1155         ENDIF
1156         !
1157         !                                                       ! non solar heat flux over the ocean (qns)
1158         IF(      srcv(jpr_qnsoce)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
1159         ELSE IF( srcv(jpr_qnsmix)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
1160         ELSE                                       ;   zqns(:,:) = 0._wp
1161         END IF
1162         ! update qns over the free ocean with:
1163         IF( nn_components /= jp_iam_opa ) THEN
1164            zqns(:,:) =  zqns(:,:) - zemp(:,:) * sst_m(:,:) * rcp         ! remove heat content due to mass flux (assumed to be at SST)
1165            IF( srcv(jpr_snow  )%laction ) THEN
1166               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean
1167            ENDIF
1168         ENDIF
1169         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
1170         ELSE                   ;   qns(:,:) =                              zqns(:,:)
1171         ENDIF
1172
1173         !                                                       ! solar flux over the ocean          (qsr)
1174         IF     ( srcv(jpr_qsroce)%laction ) THEN   ;   zqsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
1175         ELSE IF( srcv(jpr_qsrmix)%laction ) then   ;   zqsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
1176         ELSE                                       ;   zqsr(:,:) = 0._wp
1177         ENDIF
1178         IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle
1179         IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)
1180         ELSE                   ;   qsr(:,:) =                              zqsr(:,:)
1181         ENDIF
1182         !
1183         ! salt flux over the ocean (received by opa in case of opa <-> sas coupling)
1184         IF( srcv(jpr_sflx )%laction )   sfx(:,:) = frcv(jpr_sflx  )%z3(:,:,1)
1185         ! Ice cover  (received by opa in case of opa <-> sas coupling)
1186         IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
1187         !
1188
1189      ENDIF
1190      !
1191      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
1192      !
1193      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
1194      !
1195   END SUBROUTINE sbc_cpl_rcv
1196   
1197
1198   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1199      !!----------------------------------------------------------------------
1200      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
1201      !!
1202      !! ** Purpose :   provide the stress over sea-ice in coupled mode
1203      !!
1204      !! ** Method  :   transform the received stress from the atmosphere into
1205      !!             an atmosphere-ice stress in the (i,j) ocean referencial
1206      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
1207      !!                'C'-grid : i- (j-) components given at U- (V-) point
1208      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
1209      !!
1210      !!                The received stress are :
1211      !!                 - defined by 3 components (if cartesian coordinate)
1212      !!                        or by 2 components (if spherical)
1213      !!                 - oriented along geographical   coordinate (if eastward-northward)
1214      !!                        or  along the local grid coordinate (if local grid)
1215      !!                 - given at U- and V-point, resp.   if received on 2 grids
1216      !!                        or at a same point (T or I) if received on 1 grid
1217      !!                Therefore and if necessary, they are successively
1218      !!             processed in order to obtain them
1219      !!                 first  as  2 components on the sphere
1220      !!                 second as  2 components oriented along the local grid
1221      !!                 third  as  2 components on the cp_ice_msh point
1222      !!
1223      !!                Except in 'oce and ice' case, only one vector stress field
1224      !!             is received. It has already been processed in sbc_cpl_rcv
1225      !!             so that it is now defined as (i,j) components given at U-
1226      !!             and V-points, respectively. Therefore, only the third
1227      !!             transformation is done and only if the ice-grid is a 'I'-grid.
1228      !!
1229      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
1230      !!----------------------------------------------------------------------
1231      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1232      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1233      !!
1234      INTEGER ::   ji, jj                          ! dummy loop indices
1235      INTEGER ::   itx                             ! index of taux over ice
1236      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
1237      !!----------------------------------------------------------------------
1238      !
1239      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
1240      !
1241      CALL wrk_alloc( jpi,jpj, ztx, zty )
1242
1243      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
1244      ELSE                                ;   itx =  jpr_otx1
1245      ENDIF
1246
1247      ! do something only if we just received the stress from atmosphere
1248      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
1249
1250         !                                                      ! ======================= !
1251         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
1252            !                                                   ! ======================= !
1253           
1254            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
1255               !                                                       ! (cartesian to spherical -> 3 to 2 components)
1256               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
1257                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
1258               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
1259               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
1260               !
1261               IF( srcv(jpr_itx2)%laction ) THEN
1262                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
1263                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
1264                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
1265                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
1266               ENDIF
1267               !
1268            ENDIF
1269            !
1270            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
1271               !                                                       ! (geographical to local grid -> rotate the components)
1272               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
1273               IF( srcv(jpr_itx2)%laction ) THEN
1274                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
1275               ELSE
1276                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
1277               ENDIF
1278               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1279               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
1280            ENDIF
1281            !                                                   ! ======================= !
1282         ELSE                                                   !     use ocean stress    !
1283            !                                                   ! ======================= !
1284            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
1285            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
1286            !
1287         ENDIF
1288         !                                                      ! ======================= !
1289         !                                                      !     put on ice grid     !
1290         !                                                      ! ======================= !
1291         !   
1292         !                                                  j+1   j     -----V---F
1293         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
1294         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
1295         !                                                               |       |
1296         !                                                   j    j-1   -I-------|
1297         !                                               (for I)         |       |
1298         !                                                              i-1  i   i
1299         !                                                               i      i+1 (for I)
1300         SELECT CASE ( cp_ice_msh )
1301            !
1302         CASE( 'I' )                                         ! B-grid ==> I
1303            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1304            CASE( 'U' )
1305               DO jj = 2, jpjm1                                   ! (U,V) ==> I
1306                  DO ji = 2, jpim1   ! NO vector opt.
1307                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
1308                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1309                  END DO
1310               END DO
1311            CASE( 'F' )
1312               DO jj = 2, jpjm1                                   ! F ==> I
1313                  DO ji = 2, jpim1   ! NO vector opt.
1314                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
1315                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
1316                  END DO
1317               END DO
1318            CASE( 'T' )
1319               DO jj = 2, jpjm1                                   ! T ==> I
1320                  DO ji = 2, jpim1   ! NO vector opt.
1321                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
1322                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
1323                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
1324                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1325                  END DO
1326               END DO
1327            CASE( 'I' )
1328               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
1329               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1330            END SELECT
1331            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1332               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1333            ENDIF
1334            !
1335         CASE( 'F' )                                         ! B-grid ==> F
1336            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1337            CASE( 'U' )
1338               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1339                  DO ji = fs_2, fs_jpim1   ! vector opt.
1340                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1341                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1342                  END DO
1343               END DO
1344            CASE( 'I' )
1345               DO jj = 2, jpjm1                                   ! I ==> F
1346                  DO ji = 2, jpim1   ! NO vector opt.
1347                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1348                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1349                  END DO
1350               END DO
1351            CASE( 'T' )
1352               DO jj = 2, jpjm1                                   ! T ==> F
1353                  DO ji = 2, jpim1   ! NO vector opt.
1354                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1355                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1356                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1357                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1358                  END DO
1359               END DO
1360            CASE( 'F' )
1361               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1362               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1363            END SELECT
1364            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1365               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1366            ENDIF
1367            !
1368         CASE( 'C' )                                         ! C-grid ==> U,V
1369            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1370            CASE( 'U' )
1371               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1372               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1373            CASE( 'F' )
1374               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1375                  DO ji = fs_2, fs_jpim1   ! vector opt.
1376                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1377                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1378                  END DO
1379               END DO
1380            CASE( 'T' )
1381               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1382                  DO ji = fs_2, fs_jpim1   ! vector opt.
1383                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1384                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1385                  END DO
1386               END DO
1387            CASE( 'I' )
1388               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1389                  DO ji = 2, jpim1   ! NO vector opt.
1390                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1391                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1392                  END DO
1393               END DO
1394            END SELECT
1395            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1396               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1397            ENDIF
1398         END SELECT
1399
1400      ENDIF
1401      !   
1402      CALL wrk_dealloc( jpi,jpj, ztx, zty )
1403      !
1404      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
1405      !
1406   END SUBROUTINE sbc_cpl_ice_tau
1407   
1408
1409   SUBROUTINE sbc_cpl_ice_flx( p_frld, palbi, psst, pist )
1410      !!----------------------------------------------------------------------
1411      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1412      !!
1413      !! ** Purpose :   provide the heat and freshwater fluxes of the
1414      !!              ocean-ice system.
1415      !!
1416      !! ** Method  :   transform the fields received from the atmosphere into
1417      !!             surface heat and fresh water boundary condition for the
1418      !!             ice-ocean system. The following fields are provided:
1419      !!              * total non solar, solar and freshwater fluxes (qns_tot,
1420      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1421      !!             NB: emp_tot include runoffs and calving.
1422      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1423      !!             emp_ice = sublimation - solid precipitation as liquid
1424      !!             precipitation are re-routed directly to the ocean and
1425      !!             runoffs and calving directly enter the ocean.
1426      !!              * solid precipitation (sprecip), used to add to qns_tot
1427      !!             the heat lost associated to melting solid precipitation
1428      !!             over the ocean fraction.
1429      !!       ===>> CAUTION here this changes the net heat flux received from
1430      !!             the atmosphere
1431      !!
1432      !!                  - the fluxes have been separated from the stress as
1433      !!                 (a) they are updated at each ice time step compare to
1434      !!                 an update at each coupled time step for the stress, and
1435      !!                 (b) the conservative computation of the fluxes over the
1436      !!                 sea-ice area requires the knowledge of the ice fraction
1437      !!                 after the ice advection and before the ice thermodynamics,
1438      !!                 so that the stress is updated before the ice dynamics
1439      !!                 while the fluxes are updated after it.
1440      !!
1441      !! ** Action  :   update at each nf_ice time step:
1442      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1443      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1444      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
1445      !!                   emp_ice            ice sublimation - solid precipitation over the ice
1446      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
1447      !!                   sprecip             solid precipitation over the ocean 
1448      !!----------------------------------------------------------------------
1449      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1450      ! optional arguments, used only in 'mixed oce-ice' case
1451      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo
1452      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
1453      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
1454      !
1455      INTEGER ::   jl         ! dummy loop index
1456      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk
1457      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
1458      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice
1459      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3
1460      !!----------------------------------------------------------------------
1461      !
1462      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1463      !
1464      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
1465      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
1466
1467      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
1468      zicefr(:,:) = 1.- p_frld(:,:)
1469      zcptn(:,:) = rcp * sst_m(:,:)
1470      !
1471      !                                                      ! ========================= !
1472      !                                                      !    freshwater budget      !   (emp)
1473      !                                                      ! ========================= !
1474      !
1475      !                                                           ! total Precipitation - total Evaporation (emp_tot)
1476      !                                                           ! solid precipitation - sublimation       (emp_ice)
1477      !                                                           ! solid Precipitation                     (sprecip)
1478      !                                                           ! liquid + solid Precipitation            (tprecip)
1479      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1480      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1481         zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
1482         ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
1483         zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
1484         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1485            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation
1486         IF( iom_use('hflx_rain_cea') )   &
1487            CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
1488         IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   &
1489            ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1490         IF( iom_use('evap_ao_cea'  ) )   &
1491            CALL iom_put( 'evap_ao_cea'  , ztmp                   )   ! ice-free oce evap (cell average)
1492         IF( iom_use('hflx_evap_cea') )   &
1493            CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) )   ! heat flux from from evap (cell average)
1494      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1495         zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1496         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
1497         zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
1498         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
1499      END SELECT
1500
1501      IF( iom_use('subl_ai_cea') )   &
1502         CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
1503      !   
1504      !                                                           ! runoffs and calving (put in emp_tot)
1505      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1506      IF( srcv(jpr_cal)%laction ) THEN
1507         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1508         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
1509      ENDIF
1510
1511      IF( ln_mixcpl ) THEN
1512         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
1513         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
1514         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
1515         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
1516      ELSE
1517         emp_tot(:,:) =                                  zemp_tot(:,:)
1518         emp_ice(:,:) =                                  zemp_ice(:,:)
1519         sprecip(:,:) =                                  zsprecip(:,:)
1520         tprecip(:,:) =                                  ztprecip(:,:)
1521      ENDIF
1522
1523         CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
1524      IF( iom_use('snow_ao_cea') )   &
1525         CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:)             )   ! Snow        over ice-free ocean  (cell average)
1526      IF( iom_use('snow_ai_cea') )   &
1527         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
1528
1529      !                                                      ! ========================= !
1530      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1531      !                                                      ! ========================= !
1532      CASE( 'oce only' )                                     ! the required field is directly provided
1533         zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
1534      CASE( 'conservative' )                                      ! the required fields are directly provided
1535         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1536         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1537            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1538         ELSE
1539            ! Set all category values equal for the moment
1540            DO jl=1,jpl
1541               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1542            ENDDO
1543         ENDIF
1544      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
1545         zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1546         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1547            DO jl=1,jpl
1548               zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1549               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1550            ENDDO
1551         ELSE
1552            qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1553            DO jl=1,jpl
1554               zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1555               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1556            ENDDO
1557         ENDIF
1558      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
1559! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1560         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1561         zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1562            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1563            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
1564      END SELECT
1565!!gm
1566!!    currently it is taken into account in leads budget but not in the zqns_tot, and thus not in
1567!!    the flux that enter the ocean....
1568!!    moreover 1 - it is not diagnose anywhere....
1569!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1570!!
1571!! similar job should be done for snow and precipitation temperature
1572      !                                     
1573      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting
1574         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting
1575         zqns_tot(:,:) = zqns_tot(:,:) - ztmp(:,:)
1576         IF( iom_use('hflx_cal_cea') )   &
1577            CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1578      ENDIF
1579
1580      ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus
1581      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1582
1583#if defined key_lim3
1584      CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
1585
1586      ! --- evaporation --- !
1587      ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation
1588      ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice
1589      !                 but it is incoherent WITH the ice model 
1590      DO jl=1,jpl
1591         evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1)
1592      ENDDO
1593      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean
1594
1595      ! --- evaporation minus precipitation --- !
1596      emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:)
1597
1598      ! --- non solar flux over ocean --- !
1599      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
1600      zqns_oce = 0._wp
1601      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
1602
1603      ! --- heat flux associated with emp --- !
1604      zsnw(:,:) = 0._wp
1605      CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing
1606      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
1607         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip
1608         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean
1609      qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
1610         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
1611
1612      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
1613      zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
1614
1615      ! --- total non solar flux --- !
1616      zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:)
1617
1618      ! --- in case both coupled/forced are active, we must mix values --- !
1619      IF( ln_mixcpl ) THEN
1620         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
1621         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
1622         DO jl=1,jpl
1623            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
1624         ENDDO
1625         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
1626         qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
1627!!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0)
1628      ELSE
1629         qns_tot  (:,:  ) = zqns_tot  (:,:  )
1630         qns_oce  (:,:  ) = zqns_oce  (:,:  )
1631         qns_ice  (:,:,:) = zqns_ice  (:,:,:)
1632         qprec_ice(:,:)   = zqprec_ice(:,:)
1633         qemp_oce (:,:)   = zqemp_oce (:,:)
1634      ENDIF
1635
1636      CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
1637#else
1638
1639      ! clem: this formulation is certainly wrong... but better than it was...
1640      zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with:
1641         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
1642         &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
1643         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
1644
1645     IF( ln_mixcpl ) THEN
1646         qns_tot(:,:) = qns(:,:) * p_frld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
1647         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
1648         DO jl=1,jpl
1649            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
1650         ENDDO
1651      ELSE
1652         qns_tot(:,:  ) = zqns_tot(:,:  )
1653         qns_ice(:,:,:) = zqns_ice(:,:,:)
1654      ENDIF
1655
1656#endif
1657
1658      !                                                      ! ========================= !
1659      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1660      !                                                      ! ========================= !
1661      CASE( 'oce only' )
1662         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
1663      CASE( 'conservative' )
1664         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1665         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1666            zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1667         ELSE
1668            ! Set all category values equal for the moment
1669            DO jl=1,jpl
1670               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1671            ENDDO
1672         ENDIF
1673         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1674         zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1675      CASE( 'oce and ice' )
1676         zqsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1677         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1678            DO jl=1,jpl
1679               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1680               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1681            ENDDO
1682         ELSE
1683            qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1684            DO jl=1,jpl
1685               zqsr_tot(:,:   ) = zqsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1686               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1687            ENDDO
1688         ENDIF
1689      CASE( 'mixed oce-ice' )
1690         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1691! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1692!       Create solar heat flux over ice using incoming solar heat flux and albedos
1693!       ( see OASIS3 user guide, 5th edition, p39 )
1694         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1695            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1696            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1697      END SELECT
1698      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
1699         zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
1700         DO jl=1,jpl
1701            zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
1702         ENDDO
1703      ENDIF
1704
1705#if defined key_lim3
1706      CALL wrk_alloc( jpi,jpj, zqsr_oce ) 
1707      ! --- solar flux over ocean --- !
1708      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
1709      zqsr_oce = 0._wp
1710      WHERE( p_frld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / p_frld(:,:)
1711
1712      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
1713      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
1714
1715      CALL wrk_dealloc( jpi,jpj, zqsr_oce ) 
1716#endif
1717
1718      IF( ln_mixcpl ) THEN
1719         qsr_tot(:,:) = qsr(:,:) * p_frld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
1720         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
1721         DO jl=1,jpl
1722            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
1723         ENDDO
1724      ELSE
1725         qsr_tot(:,:  ) = zqsr_tot(:,:  )
1726         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
1727      ENDIF
1728
1729      !                                                      ! ========================= !
1730      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
1731      !                                                      ! ========================= !
1732      CASE ('coupled')
1733         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1734            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
1735         ELSE
1736            ! Set all category values equal for the moment
1737            DO jl=1,jpl
1738               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
1739            ENDDO
1740         ENDIF
1741      END SELECT
1742     
1743      IF( ln_mixcpl ) THEN
1744         DO jl=1,jpl
1745            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
1746         ENDDO
1747      ELSE
1748         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
1749      ENDIF
1750     
1751      !                                                      ! ========================= !
1752      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
1753      !                                                      ! ========================= !
1754      CASE ('coupled')
1755         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
1756         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
1757      END SELECT
1758
1759      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
1760      ! Used for LIM2 and LIM3
1761      ! Coupled case: since cloud cover is not received from atmosphere
1762      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
1763      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
1764      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
1765
1766      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
1767      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
1768      !
1769      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
1770      !
1771   END SUBROUTINE sbc_cpl_ice_flx
1772   
1773   
1774   SUBROUTINE sbc_cpl_snd( kt )
1775      !!----------------------------------------------------------------------
1776      !!             ***  ROUTINE sbc_cpl_snd  ***
1777      !!
1778      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
1779      !!
1780      !! ** Method  :   send to the atmosphere through a call to cpl_snd
1781      !!              all the needed fields (as defined in sbc_cpl_init)
1782      !!----------------------------------------------------------------------
1783      INTEGER, INTENT(in) ::   kt
1784      !
1785      INTEGER ::   ji, jj, jl   ! dummy loop indices
1786      INTEGER ::   isec, info   ! local integer
1787      REAL(wp) ::   zumax, zvmax
1788      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
1789      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
1790      !!----------------------------------------------------------------------
1791      !
1792      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
1793      !
1794      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1795      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1796
1797      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
1798
1799      zfr_l(:,:) = 1.- fr_i(:,:)
1800      !                                                      ! ------------------------- !
1801      !                                                      !    Surface temperature    !   in Kelvin
1802      !                                                      ! ------------------------- !
1803      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
1804         
1805         IF ( nn_components == jp_iam_opa ) THEN
1806            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of ln_useCT on the received part
1807         ELSE
1808            ! we must send the surface potential temperature
1809            IF( ln_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
1810            ELSE                    ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
1811            ENDIF
1812            !
1813            SELECT CASE( sn_snd_temp%cldes)
1814            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
1815            CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
1816               SELECT CASE( sn_snd_temp%clcat )
1817               CASE( 'yes' )   
1818                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl)
1819               CASE( 'no' )
1820                  WHERE( SUM( a_i, dim=3 ) /= 0. )
1821                     ztmp3(:,:,1) = SUM( tn_ice * a_i, dim=3 ) / SUM( a_i, dim=3 )
1822                  ELSEWHERE
1823                     ztmp3(:,:,1) = rt0 ! TODO: Is freezing point a good default? (Maybe SST is better?)
1824                  END WHERE
1825               CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1826               END SELECT
1827            CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)   
1828               SELECT CASE( sn_snd_temp%clcat )
1829               CASE( 'yes' )   
1830                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1831               CASE( 'no' )
1832                  ztmp3(:,:,:) = 0.0
1833                  DO jl=1,jpl
1834                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
1835                  ENDDO
1836               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1837               END SELECT
1838            CASE( 'mixed oce-ice'        )   
1839               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) 
1840               DO jl=1,jpl
1841                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
1842               ENDDO
1843            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
1844            END SELECT
1845         ENDIF
1846         IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1847         IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
1848         IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1849      ENDIF
1850      !                                                      ! ------------------------- !
1851      !                                                      !           Albedo          !
1852      !                                                      ! ------------------------- !
1853      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
1854         SELECT CASE( sn_snd_alb%cldes )
1855         CASE( 'ice'          )   ; ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
1856         CASE( 'weighted ice' )   ; ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1857         CASE default             ; CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
1858         END SELECT
1859         CALL cpl_snd( jps_albice, isec, ztmp3, info )
1860      ENDIF
1861      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
1862         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
1863         DO jl=1,jpl
1864            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
1865         ENDDO
1866         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1867      ENDIF
1868      !                                                      ! ------------------------- !
1869      !                                                      !  Ice fraction & Thickness !
1870      !                                                      ! ------------------------- !
1871      ! Send ice fraction field to atmosphere
1872      IF( ssnd(jps_fice)%laction ) THEN
1873         SELECT CASE( sn_snd_thick%clcat )
1874         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
1875         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
1876         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1877         END SELECT
1878         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
1879      ENDIF
1880     
1881      ! Send ice fraction field to OPA (sent by SAS in SAS-OPA coupling)
1882      IF( ssnd(jps_fice2)%laction ) THEN
1883         ztmp3(:,:,1) = fr_i(:,:)
1884         IF( ssnd(jps_fice2)%laction )   CALL cpl_snd( jps_fice2, isec, ztmp3, info )
1885      ENDIF
1886
1887      ! Send ice and snow thickness field
1888      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
1889         SELECT CASE( sn_snd_thick%cldes)
1890         CASE( 'none'                  )       ! nothing to do
1891         CASE( 'weighted ice and snow' )   
1892            SELECT CASE( sn_snd_thick%clcat )
1893            CASE( 'yes' )   
1894               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
1895               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
1896            CASE( 'no' )
1897               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
1898               DO jl=1,jpl
1899                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
1900                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
1901               ENDDO
1902            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1903            END SELECT
1904         CASE( 'ice and snow'         )   
1905            SELECT CASE( sn_snd_thick%clcat )
1906            CASE( 'yes' )
1907               ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
1908               ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
1909            CASE( 'no' )
1910               WHERE( SUM( a_i, dim=3 ) /= 0. )
1911                  ztmp3(:,:,1) = SUM( ht_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
1912                  ztmp4(:,:,1) = SUM( ht_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
1913               ELSEWHERE
1914                 ztmp3(:,:,1) = 0.
1915                 ztmp4(:,:,1) = 0.
1916               END WHERE
1917            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1918            END SELECT
1919         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
1920         END SELECT
1921         IF( ssnd(jps_hice)%laction )   CALL cpl_snd( jps_hice, isec, ztmp3, info )
1922         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
1923      ENDIF
1924      !
1925#if defined key_cpl_carbon_cycle
1926      !                                                      ! ------------------------- !
1927      !                                                      !  CO2 flux from PISCES     !
1928      !                                                      ! ------------------------- !
1929      IF( ssnd(jps_co2)%laction )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
1930      !
1931#endif
1932      !                                                      ! ------------------------- !
1933      !                                                      !      Surface current      !
1934      !                                                      ! ------------------------- !
1935      IF( ssnd(jps_ocx1)%laction ) THEN
1936         !   
1937         !                                                  j+1   j     -----V---F
1938         ! surface velocity always sent from T point                     !       |
1939         !                                                        j      |   T   U
1940         !                                                               |       |
1941         !                                                   j    j-1   -I-------|
1942         !                                               (for I)         |       |
1943         !                                                              i-1  i   i
1944         !                                                               i      i+1 (for I)
1945         IF( nn_components == jp_iam_opa ) THEN
1946            zotx1(:,:) = un(:,:,1) 
1947            zoty1(:,:) = vn(:,:,1) 
1948         ELSE       
1949            SELECT CASE( TRIM( sn_snd_crt%cldes ) )
1950            CASE( 'oce only'             )      ! C-grid ==> T
1951               DO jj = 2, jpjm1
1952                  DO ji = fs_2, fs_jpim1   ! vector opt.
1953                     zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
1954                     zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
1955                  END DO
1956               END DO
1957            CASE( 'weighted oce and ice' )   
1958               SELECT CASE ( cp_ice_msh )
1959               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1960                  DO jj = 2, jpjm1
1961                     DO ji = fs_2, fs_jpim1   ! vector opt.
1962                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1963                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
1964                        zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1965                        zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1966                     END DO
1967                  END DO
1968               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1969                  DO jj = 2, jpjm1
1970                     DO ji = 2, jpim1   ! NO vector opt.
1971                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1972                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1973                        zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1974                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1975                        zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1976                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1977                     END DO
1978                  END DO
1979               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1980                  DO jj = 2, jpjm1
1981                     DO ji = 2, jpim1   ! NO vector opt.
1982                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1983                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1984                        zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1985                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1986                        zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1987                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1988                     END DO
1989                  END DO
1990               END SELECT
1991               CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
1992            CASE( 'mixed oce-ice'        )
1993               SELECT CASE ( cp_ice_msh )
1994               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1995                  DO jj = 2, jpjm1
1996                     DO ji = fs_2, fs_jpim1   ! vector opt.
1997                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
1998                           &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1999                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
2000                           &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2001                     END DO
2002                  END DO
2003               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2004                  DO jj = 2, jpjm1
2005                     DO ji = 2, jpim1   ! NO vector opt.
2006                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2007                           &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2008                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2009                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2010                           &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2011                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2012                     END DO
2013                  END DO
2014               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2015                  DO jj = 2, jpjm1
2016                     DO ji = 2, jpim1   ! NO vector opt.
2017                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2018                           &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2019                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2020                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2021                           &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2022                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2023                     END DO
2024                  END DO
2025               END SELECT
2026            END SELECT
2027            CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
2028            !
2029         ENDIF
2030         !
2031         !
2032         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2033            !                                                                     ! Ocean component
2034            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
2035            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
2036            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
2037            zoty1(:,:) = ztmp2(:,:)
2038            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
2039               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
2040               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
2041               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
2042               zity1(:,:) = ztmp2(:,:)
2043            ENDIF
2044         ENDIF
2045         !
2046         ! spherical coordinates to cartesian -> 2 components to 3 components
2047         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
2048            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2049            ztmp2(:,:) = zoty1(:,:)
2050            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2051            !
2052            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2053               ztmp1(:,:) = zitx1(:,:)
2054               ztmp1(:,:) = zity1(:,:)
2055               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2056            ENDIF
2057         ENDIF
2058         !
2059         IF( ssnd(jps_ocx1)%laction )   CALL cpl_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2060         IF( ssnd(jps_ocy1)%laction )   CALL cpl_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2061         IF( ssnd(jps_ocz1)%laction )   CALL cpl_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
2062         !
2063         IF( ssnd(jps_ivx1)%laction )   CALL cpl_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
2064         IF( ssnd(jps_ivy1)%laction )   CALL cpl_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
2065         IF( ssnd(jps_ivz1)%laction )   CALL cpl_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
2066      !
2067      ENDIF
2068      !                                                      ! ------------------------- !
2069      !                                                      !  Surface current to waves !
2070      !                                                      ! ------------------------- !
2071      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN
2072         !   
2073         !                                                  j+1   j     -----V---F
2074         ! surface velocity always sent from T point                     !       |
2075         !                                                        j      |   T   U
2076         !                                                               |       |
2077         !                                                   j    j-1   -I-------|
2078         !                                               (for I)         |       |
2079         !                                                              i-1  i   i
2080         !                                                               i      i+1 (for I)
2081         SELECT CASE( TRIM( sn_snd_crtw%cldes ) )
2082         CASE( 'oce only'             )      ! C-grid ==> T
2083            DO jj = 2, jpjm1
2084               DO ji = fs_2, fs_jpim1   ! vector opt.
2085                  zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
2086                  zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
2087               END DO
2088            END DO
2089         CASE( 'weighted oce and ice' )   
2090            SELECT CASE ( cp_ice_msh )
2091            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2092               DO jj = 2, jpjm1
2093                  DO ji = fs_2, fs_jpim1   ! vector opt.
2094                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2095                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
2096                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2097                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2098                  END DO
2099               END DO
2100            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2101               DO jj = 2, jpjm1
2102                  DO ji = 2, jpim1   ! NO vector opt.
2103                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2104                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2105                     zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2106                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2107                     zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2108                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2109                  END DO
2110               END DO
2111            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2112               DO jj = 2, jpjm1
2113                  DO ji = 2, jpim1   ! NO vector opt.
2114                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2115                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2116                     zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2117                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2118                     zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2119                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2120                  END DO
2121               END DO
2122            END SELECT
2123            CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
2124         CASE( 'mixed oce-ice'        )
2125            SELECT CASE ( cp_ice_msh )
2126            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2127               DO jj = 2, jpjm1
2128                  DO ji = fs_2, fs_jpim1   ! vector opt.
2129                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
2130                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2131                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
2132                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2133                  END DO
2134               END DO
2135            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2136               DO jj = 2, jpjm1
2137                  DO ji = 2, jpim1   ! NO vector opt.
2138                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2139                        &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2140                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2141                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2142                        &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2143                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2144                  END DO
2145               END DO
2146            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2147               DO jj = 2, jpjm1
2148                  DO ji = 2, jpim1   ! NO vector opt.
2149                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2150                        &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2151                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2152                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2153                        &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2154                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2155                  END DO
2156               END DO
2157            END SELECT
2158         END SELECT
2159         CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. )
2160         !
2161         !
2162         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2163            !                                                                     ! Ocean component
2164            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component
2165            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
2166            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
2167            zoty1(:,:) = ztmp2(:,:)
2168            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
2169               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
2170               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
2171               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
2172               zity1(:,:) = ztmp2(:,:)
2173            ENDIF
2174         ENDIF
2175         !
2176!         ! spherical coordinates to cartesian -> 2 components to 3 components
2177!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
2178!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2179!            ztmp2(:,:) = zoty1(:,:)
2180!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2181!            !
2182!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2183!               ztmp1(:,:) = zitx1(:,:)
2184!               ztmp1(:,:) = zity1(:,:)
2185!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2186!            ENDIF
2187!         ENDIF
2188         !
2189         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2190         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2191      !
2192      ENDIF
2193      !
2194      IF( ssnd(jps_ficet)%laction ) THEN
2195           CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info )
2196      END IF
2197      !                                                      ! ------------------------- !
2198      !                                                      !   Water levels to waves   !
2199      !                                                      ! ------------------------- !
2200      IF( ssnd(jps_wlev)%laction ) THEN
2201           IF( ln_apr_dyn ) THEN
2202               IF( kt /= nit000 ) THEN
2203                   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
2204               ELSE
2205                   ztmp1(:,:) = sshb(:,:) 
2206               ENDIF
2207           ELSE
2208               ztmp1(:,:) = sshn(:,:) 
2209           ENDIF 
2210           ! add the bathymetry - we need to send the total height of the water column
2211           ztmp1(:,:) = ztmp1(:,:) + bathy(:,:) 
2212           CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2213      END IF
2214      !
2215      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
2216      !                                                        ! SSH
2217      IF( ssnd(jps_ssh )%laction )  THEN
2218         !                          ! removed inverse barometer ssh when Patm
2219         !                          forcing is used (for sea-ice dynamics)
2220         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
2221         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
2222         ENDIF
2223         CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
2224      ENDIF
2225      !                                                        ! SSS
2226      IF( ssnd(jps_soce  )%laction )  THEN
2227         CALL cpl_snd( jps_soce  , isec, RESHAPE ( tsn(:,:,1,jp_sal), (/jpi,jpj,1/) ), info )
2228      ENDIF
2229      !                                                        ! first T level thickness
2230      IF( ssnd(jps_e3t1st )%laction )  THEN
2231         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( fse3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
2232      ENDIF
2233      !                                                        ! Qsr fraction
2234      IF( ssnd(jps_fraqsr)%laction )  THEN
2235         CALL cpl_snd( jps_fraqsr, isec, RESHAPE ( fraqsr_1lev(:,:) , (/jpi,jpj,1/) ), info )
2236      ENDIF
2237      !
2238      !  Fields sent by SAS to OPA when OASIS coupling
2239      !                                                        ! Solar heat flux
2240      IF( ssnd(jps_qsroce)%laction )  CALL cpl_snd( jps_qsroce, isec, RESHAPE ( qsr , (/jpi,jpj,1/) ), info )
2241      IF( ssnd(jps_qnsoce)%laction )  CALL cpl_snd( jps_qnsoce, isec, RESHAPE ( qns , (/jpi,jpj,1/) ), info )
2242      IF( ssnd(jps_oemp  )%laction )  CALL cpl_snd( jps_oemp  , isec, RESHAPE ( emp , (/jpi,jpj,1/) ), info )
2243      IF( ssnd(jps_sflx  )%laction )  CALL cpl_snd( jps_sflx  , isec, RESHAPE ( sfx , (/jpi,jpj,1/) ), info )
2244      IF( ssnd(jps_otx1  )%laction )  CALL cpl_snd( jps_otx1  , isec, RESHAPE ( utau, (/jpi,jpj,1/) ), info )
2245      IF( ssnd(jps_oty1  )%laction )  CALL cpl_snd( jps_oty1  , isec, RESHAPE ( vtau, (/jpi,jpj,1/) ), info )
2246      IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info )
2247      IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
2248
2249      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
2250      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
2251      !
2252      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
2253      !
2254   END SUBROUTINE sbc_cpl_snd
2255   
2256   !!======================================================================
2257END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.