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

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

source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 3875

Last change on this file since 3875 was 3875, checked in by clevy, 11 years ago

Configuration Setting/Step? 1, see ticket:#1074

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