source: branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5279

Last change on this file since 5279 was 5279, checked in by dancopsey, 6 years ago

Removed SVN keywords.

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