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

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

source: branches/2012/dev_LOCEAN_2012/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 3647

Last change on this file since 3647 was 3647, checked in by smasson, 11 years ago

coupling interface without sea-ice, see ticket:1018

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