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

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

source: branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5490

Last change on this file since 5490 was 5490, checked in by dancopsey, 9 years ago

Reintroduced acidentally deleted CASE command.

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