source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils305/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 12003

Last change on this file since 12003 was 12003, checked in by dford, 10 months ago

Merge in change relating to pressure correction.

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