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/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5236

Last change on this file since 5236 was 5236, checked in by cetlod, 9 years ago

NEMOGCM_dev_r5204_CNRS_PISCES_dcy : update routines according to the new strategy, see ticket #1484

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