source: branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5222

Last change on this file since 5222 was 5222, checked in by cetlod, 6 years ago

dev_r5204_CNRS_PISCES_dcy: some improvments

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