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

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

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

Last change on this file since 4161 was 4161, checked in by cetlod, 9 years ago

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

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