source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 3402

Last change on this file since 3402 was 3402, checked in by acc, 9 years ago

Branch: dev_r3385_NOCS04_HAMF; #665. Stage 2 of 2012 development: suppression of emps array and introduction of sfx (salt flux) array with associated code to setup the options for embedding the seaice into the ocean

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