source: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbccpl.F90 @ 2194

Last change on this file since 2194 was 2194, checked in by aclsce, 10 years ago

Modified NEMO/SOURCES specific to IPSLCM6 configuration ( compilation files and NEMO model sources)..
Removed some routines which are not specific anymore to IPSL configuration.

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