source: branches/2013/dev_r4022_MERCATOR5_PAPA1D/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 4210

Last change on this file since 4210 was 4210, checked in by cbricaud, 8 years ago

merge changes from rev 4022 to 4119 from trunk into dev_r4022_MERCATOR5_PAPA1D

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