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

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

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 4069

Last change on this file since 4069 was 4069, checked in by charris, 11 years ago

#1107 Change as on trunk.

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