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

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

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5033

Last change on this file since 5033 was 5030, checked in by timgraham, 9 years ago

Added changes for coupling of meltpond fraction and depth between NEMO-CICE and atmosphere model

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