source: LMDZ5/trunk/libf/phylmd/cpl_mod.F90 @ 2841

Last change on this file since 2841 was 2545, checked in by Laurent Fairhead, 8 years ago

Rollback for revision r2538: calculation of the fluxes is correct but transmission
to coupler is broken
LF

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 55.9 KB
Line 
1!
2MODULE cpl_mod
3!
4! This module excahanges and transforms all fields that should be recieved or sent to
5! coupler. The transformation of the fields are done from the grid 1D-array in phylmd
6! to the regular 2D grid accepted by the coupler. Cumulation of the fields for each
7! timestep is done in here.
8!
9! Each type of surface that recevie fields from the coupler have a subroutine named
10! cpl_receive_XXX_fields and each surface that have fields to be sent to the coupler
11! have a subroutine named cpl_send_XXX_fields.
12!
13!*************************************************************************************
14
15! Use statements
16!*************************************************************************************
17  USE dimphy, ONLY : klon
18  USE mod_phys_lmdz_para
19  USE ioipsl
20  USE iophy
21
22! The module oasis is always used. Without the cpp key CPP_COUPLE only the parameters
23! in the module are compiled and not the subroutines.
24  USE oasis
25  USE write_field_phy
26  USE time_phylmdz_mod, ONLY: day_step_phy
27 
28! Global attributes
29!*************************************************************************************
30  IMPLICIT NONE
31  PRIVATE
32
33  ! All subroutine are public except cpl_send_all
34  PUBLIC :: cpl_init, cpl_receive_frac, cpl_receive_ocean_fields, cpl_receive_seaice_fields, &
35       cpl_send_ocean_fields, cpl_send_seaice_fields, cpl_send_land_fields, &
36       cpl_send_landice_fields, gath2cpl
37 
38
39! Declaration of module variables
40!*************************************************************************************
41! variable for coupling period
42  INTEGER, SAVE :: nexca
43  !$OMP THREADPRIVATE(nexca)
44
45! variables for cumulating fields during a coupling periode :
46  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_sols, cpl_nsol, cpl_rain
47  !$OMP THREADPRIVATE(cpl_sols,cpl_nsol,cpl_rain)
48  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_snow, cpl_evap, cpl_tsol
49  !$OMP THREADPRIVATE(cpl_snow,cpl_evap,cpl_tsol)
50  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_fder, cpl_albe, cpl_taux, cpl_tauy
51  !$OMP THREADPRIVATE(cpl_fder,cpl_albe,cpl_taux,cpl_tauy)
52  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_windsp
53  !$OMP THREADPRIVATE(cpl_windsp)
54  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_taumod
55  !$OMP THREADPRIVATE(cpl_taumod)
56  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_atm_co2
57  !$OMP THREADPRIVATE(cpl_atm_co2)
58  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_rriv2D, cpl_rcoa2D, cpl_rlic2D
59  !$OMP THREADPRIVATE(cpl_rriv2D,cpl_rcoa2D,cpl_rlic2D)
60
61! variables read from coupler :
62  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_sst     ! sea surface temperature
63  !$OMP THREADPRIVATE(read_sst)
64  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_sit     ! sea ice temperature
65  !$OMP THREADPRIVATE(read_sit)
66  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_sic     ! sea ice fraction
67  !$OMP THREADPRIVATE(read_sic)
68  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_alb_sic ! albedo at sea ice
69  !$OMP THREADPRIVATE(read_alb_sic)
70  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_u0, read_v0 ! ocean surface current
71  !$OMP THREADPRIVATE(read_u0,read_v0)
72  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_co2     ! ocean co2 flux
73  !$OMP THREADPRIVATE(read_co2)
74  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: unity
75  !$OMP THREADPRIVATE(unity)
76  INTEGER, SAVE                             :: nidct, nidcs
77  !$OMP THREADPRIVATE(nidct,nidcs)
78
79! variables to be sent to the coupler
80  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_sols2D, cpl_nsol2D, cpl_rain2D
81  !$OMP THREADPRIVATE(cpl_sols2D, cpl_nsol2D, cpl_rain2D)
82  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_snow2D, cpl_evap2D, cpl_tsol2D
83  !$OMP THREADPRIVATE(cpl_snow2D, cpl_evap2D, cpl_tsol2D)
84  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_fder2D, cpl_albe2D
85  !$OMP THREADPRIVATE(cpl_fder2D, cpl_albe2D)
86  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_taux2D, cpl_tauy2D
87  !$OMP THREADPRIVATE(cpl_taux2D, cpl_tauy2D)
88  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_taumod2D
89  !$OMP THREADPRIVATE(cpl_taumod2D)
90  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_windsp2D
91  !$OMP THREADPRIVATE(cpl_windsp2D)
92  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_atm_co22D
93  !$OMP THREADPRIVATE(cpl_atm_co22D)
94
95CONTAINS
96!
97!************************************************************************************
98!
99  SUBROUTINE cpl_init(dtime, rlon, rlat)
100    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_ocn_day
101    USE surface_data
102    USE indice_sol_mod
103    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
104    USE time_phylmdz_mod, ONLY: annee_ref, day_ini, itau_phy, itaufin_phy
105    USE print_control_mod, ONLY: lunout
106
107! Input arguments
108!*************************************************************************************
109    REAL, INTENT(IN)                  :: dtime
110    REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
111
112! Local variables
113!*************************************************************************************
114    INTEGER                           :: error, sum_error, ig, i
115    INTEGER                           :: jf, nhoridct
116    INTEGER                           :: nhoridcs
117    INTEGER                           :: idtime
118    INTEGER                           :: idayref
119    INTEGER                           :: npas ! only for OASIS2
120    REAL                              :: zjulian
121    REAL, DIMENSION(nbp_lon,nbp_lat)  :: zx_lon, zx_lat
122    CHARACTER(len = 20)               :: modname = 'cpl_init'
123    CHARACTER(len = 80)               :: abort_message
124    CHARACTER(len=80)                 :: clintocplnam, clfromcplnam
125
126!*************************************************************************************
127! Calculate coupling period
128!
129!*************************************************************************************
130     
131    npas = itaufin_phy
132!    nexca = 86400 / dtime
133    nexca = t_coupl / dtime
134    WRITE(lunout,*)' ##### Ocean couple #####'
135    WRITE(lunout,*)' Valeurs des pas de temps'
136    WRITE(lunout,*)' npas = ', npas
137    WRITE(lunout,*)' nexca = ', nexca
138   
139!*************************************************************************************
140! Allocate variables
141!
142!*************************************************************************************
143    error = 0
144    sum_error = 0
145
146    ALLOCATE(unity(klon), stat = error)
147    sum_error = sum_error + error
148    ALLOCATE(cpl_sols(klon,2), stat = error)
149    sum_error = sum_error + error
150    ALLOCATE(cpl_nsol(klon,2), stat = error)
151    sum_error = sum_error + error
152    ALLOCATE(cpl_rain(klon,2), stat = error)
153    sum_error = sum_error + error
154    ALLOCATE(cpl_snow(klon,2), stat = error)
155    sum_error = sum_error + error
156    ALLOCATE(cpl_evap(klon,2), stat = error)
157    sum_error = sum_error + error
158    ALLOCATE(cpl_tsol(klon,2), stat = error)
159    sum_error = sum_error + error
160    ALLOCATE(cpl_fder(klon,2), stat = error)
161    sum_error = sum_error + error
162    ALLOCATE(cpl_albe(klon,2), stat = error)
163    sum_error = sum_error + error
164    ALLOCATE(cpl_taux(klon,2), stat = error)
165    sum_error = sum_error + error
166    ALLOCATE(cpl_tauy(klon,2), stat = error)
167    sum_error = sum_error + error
168    ALLOCATE(cpl_windsp(klon,2), stat = error)
169    sum_error = sum_error + error
170    ALLOCATE(cpl_taumod(klon,2), stat = error)
171    sum_error = sum_error + error
172    ALLOCATE(cpl_rriv2D(nbp_lon,jj_nb), stat=error)
173    sum_error = sum_error + error
174    ALLOCATE(cpl_rcoa2D(nbp_lon,jj_nb), stat=error)
175    sum_error = sum_error + error
176    ALLOCATE(cpl_rlic2D(nbp_lon,jj_nb), stat=error)
177    sum_error = sum_error + error
178    ALLOCATE(read_sst(nbp_lon, jj_nb), stat = error)
179    sum_error = sum_error + error
180    ALLOCATE(read_sic(nbp_lon, jj_nb), stat = error)
181    sum_error = sum_error + error
182    ALLOCATE(read_sit(nbp_lon, jj_nb), stat = error)
183    sum_error = sum_error + error
184    ALLOCATE(read_alb_sic(nbp_lon, jj_nb), stat = error)
185    sum_error = sum_error + error
186    ALLOCATE(read_u0(nbp_lon, jj_nb), stat = error)
187    sum_error = sum_error + error
188    ALLOCATE(read_v0(nbp_lon, jj_nb), stat = error)
189    sum_error = sum_error + error
190
191    IF (carbon_cycle_cpl) THEN
192       ALLOCATE(read_co2(nbp_lon, jj_nb), stat = error)
193       sum_error = sum_error + error
194       ALLOCATE(cpl_atm_co2(klon,2), stat = error)
195       sum_error = sum_error + error
196
197! Allocate variable in carbon_cycle_mod
198       ALLOCATE(fco2_ocn_day(klon), stat = error)
199       sum_error = sum_error + error
200    END IF
201
202    IF (sum_error /= 0) THEN
203       abort_message='Pb allocation variables couplees'
204       CALL abort_physic(modname,abort_message,1)
205    ENDIF
206!*************************************************************************************
207! Initialize the allocated varaibles
208!
209!*************************************************************************************
210    DO ig = 1, klon
211       unity(ig) = ig
212    ENDDO
213
214!*************************************************************************************
215! Initialize coupling
216!
217!*************************************************************************************
218    idtime = INT(dtime)
219#ifdef CPP_COUPLE
220    CALL inicma
221#endif
222
223!*************************************************************************************
224! initialize NetCDF output
225!
226!*************************************************************************************
227    IF (is_sequential) THEN
228       idayref = day_ini
229       CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
230       CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlon,zx_lon)
231       DO i = 1, nbp_lon
232          zx_lon(i,1) = rlon(i+1)
233          zx_lon(i,nbp_lat) = rlon(i+1)
234       ENDDO
235       CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlat,zx_lat)
236       clintocplnam="cpl_atm_tauflx"
237       CALL histbeg(clintocplnam,nbp_lon,zx_lon(:,1),nbp_lat,zx_lat(1,:),&
238            1,nbp_lon,1,nbp_lat, itau_phy,zjulian,dtime,nhoridct,nidct)
239! no vertical axis
240       CALL histdef(nidct, 'tauxe','tauxe', &
241            "-",nbp_lon,nbp_lat, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
242       CALL histdef(nidct, 'tauyn','tauyn', &
243            "-",nbp_lon,nbp_lat, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
244       CALL histdef(nidct, 'tmp_lon','tmp_lon', &
245            "-",nbp_lon,nbp_lat, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
246       CALL histdef(nidct, 'tmp_lat','tmp_lat', &
247            "-",nbp_lon,nbp_lat, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
248       DO jf=1,maxsend
249         IF (infosend(i)%action) THEN
250             CALL histdef(nidct, infosend(i)%name ,infosend(i)%name , &
251                "-",nbp_lon,nbp_lat,nhoridct,1,1,1,-99,32,"inst",dtime,dtime)
252         ENDIF
253       END DO
254       CALL histend(nidct)
255       CALL histsync(nidct)
256       
257       clfromcplnam="cpl_atm_sst"
258       CALL histbeg(clfromcplnam,nbp_lon,zx_lon(:,1),nbp_lat,zx_lat(1,:),1,nbp_lon,1,nbp_lat, &
259            0,zjulian,dtime,nhoridcs,nidcs)
260! no vertical axis
261       DO jf=1,maxrecv
262         IF (inforecv(i)%action) THEN
263             CALL histdef(nidcs,inforecv(i)%name ,inforecv(i)%name , &
264                "-",nbp_lon,nbp_lat,nhoridcs,1,1,1,-99,32,"inst",dtime,dtime)
265         ENDIF
266       END DO
267       CALL histend(nidcs)
268       CALL histsync(nidcs)
269
270    ENDIF    ! is_sequential
271   
272
273!*************************************************************************************
274! compatibility test
275!
276!*************************************************************************************
277    IF (carbon_cycle_cpl .AND. version_ocean=='opa8') THEN
278       abort_message='carbon_cycle_cpl does not work with opa8'
279       CALL abort_physic(modname,abort_message,1)
280    END IF
281
282  END SUBROUTINE cpl_init
283 
284!
285!*************************************************************************************
286!
287 
288  SUBROUTINE cpl_receive_frac(itime, dtime, pctsrf, is_modified)
289! This subroutine receives from coupler for both ocean and seaice
290! 4 fields : read_sst, read_sic, read_sit and read_alb_sic.
291! The new sea-ice-land-landice fraction is returned. The others fields
292! are stored in this module.
293    USE surface_data
294    USE geometry_mod, ONLY : longitude_deg, latitude_deg
295    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl
296    USE indice_sol_mod
297    USE time_phylmdz_mod, ONLY: start_time, itau_phy
298    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
299
300    INCLUDE "YOMCST.h"
301
302! Arguments
303!************************************************************************************
304    INTEGER, INTENT(IN)                        :: itime
305    REAL, INTENT(IN)                           :: dtime
306    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf
307    LOGICAL, INTENT(OUT)                       :: is_modified
308
309! Local variables
310!************************************************************************************
311    INTEGER                                 :: j, i, time_sec
312    INTEGER                                 :: itau_w
313    INTEGER, DIMENSION(nbp_lon*nbp_lat)     :: ndexcs
314    CHARACTER(len = 20)                     :: modname = 'cpl_receive_frac'
315    CHARACTER(len = 80)                     :: abort_message
316    REAL, DIMENSION(klon)                   :: read_sic1D
317    REAL, DIMENSION(nbp_lon,jj_nb,maxrecv)      :: tab_read_flds
318    REAL, DIMENSION(klon,nbsrf)             :: pctsrf_old
319    REAL, DIMENSION(klon_mpi)               :: rlon_mpi, rlat_mpi
320    REAL, DIMENSION(nbp_lon,jj_nb)             :: tmp_lon, tmp_lat
321    REAL, DIMENSION(nbp_lon,jj_nb)             :: tmp_r0
322
323!*************************************************************************************
324! Start calculation
325! Get fields from coupler
326!
327!*************************************************************************************
328
329    is_modified=.FALSE.
330
331! Check if right moment to receive from coupler
332    IF (MOD(itime, nexca) == 1) THEN
333       is_modified=.TRUE.
334 
335       time_sec=(itime-1)*dtime
336#ifdef CPP_COUPLE
337!$OMP MASTER
338    CALL fromcpl(time_sec, tab_read_flds)
339!$OMP END MASTER
340#endif
341   
342! NetCDF output of received fields
343       IF (is_sequential) THEN
344          ndexcs(:) = 0
345          itau_w = itau_phy + itime + start_time * day_step_phy
346          DO i = 1, maxrecv
347            IF (inforecv(i)%action) THEN
348                CALL histwrite(nidcs,inforecv(i)%name,itau_w,tab_read_flds(:,:,i),nbp_lon*(nbp_lat),ndexcs)
349            ENDIF
350          END DO
351       ENDIF
352
353
354! Save each field in a 2D array.
355!$OMP MASTER
356       read_sst(:,:)     = tab_read_flds(:,:,idr_sisutw)  ! Sea surface temperature
357       read_sic(:,:)     = tab_read_flds(:,:,idr_icecov)  ! Sea ice concentration
358       read_alb_sic(:,:) = tab_read_flds(:,:,idr_icealw)  ! Albedo at sea ice
359       read_sit(:,:)     = tab_read_flds(:,:,idr_icetem)  ! Sea ice temperature
360!$OMP END MASTER
361
362       IF (cpl_current) THEN
363
364! Transform the longitudes and latitudes on 2D arrays
365          CALL gather_omp(longitude_deg,rlon_mpi)
366          CALL gather_omp(latitude_deg,rlat_mpi)
367!$OMP MASTER
368          CALL Grid1DTo2D_mpi(rlon_mpi,tmp_lon)
369          CALL Grid1DTo2D_mpi(rlat_mpi,tmp_lat)
370
371! Transform the currents from cartesian to spheric coordinates
372! tmp_r0 should be zero
373          CALL geo2atm(nbp_lon, jj_nb, tab_read_flds(:,:,idr_curenx), &
374             tab_read_flds(:,:,idr_cureny), tab_read_flds(:,:,idr_curenz), &
375               tmp_lon, tmp_lat, &
376               read_u0(:,:), read_v0(:,:), tmp_r0(:,:))
377!$OMP END MASTER
378
379      ELSE
380          read_u0(:,:) = 0.
381          read_v0(:,:) = 0.
382      ENDIF
383
384       IF (carbon_cycle_cpl) THEN
385!$OMP MASTER
386           read_co2(:,:) = tab_read_flds(:,:,idr_oceco2) ! CO2 flux
387!$OMP END MASTER
388       ENDIF
389
390!*************************************************************************************
391!  Transform seaice fraction (read_sic : ocean-seaice mask) into global
392!  fraction (pctsrf : ocean-seaice-land-landice mask)
393!
394!*************************************************************************************
395       CALL cpl2gath(read_sic, read_sic1D, klon, unity)
396
397       pctsrf_old(:,:) = pctsrf(:,:)
398       DO i = 1, klon
399          ! treatment only of points with ocean and/or seaice
400          ! old land-ocean mask can not be changed
401          IF (pctsrf_old(i,is_oce) + pctsrf_old(i,is_sic) > 0.) THEN
402             pctsrf(i,is_sic) = (pctsrf_old(i,is_oce) + pctsrf_old(i,is_sic)) &
403                  * read_sic1D(i)
404             pctsrf(i,is_oce) = (pctsrf_old(i,is_oce) + pctsrf_old(i,is_sic)) &
405                  - pctsrf(i,is_sic)
406          ENDIF
407       ENDDO
408
409    END IF ! if time to receive
410
411  END SUBROUTINE cpl_receive_frac
412
413!
414!*************************************************************************************
415!
416
417  SUBROUTINE cpl_receive_ocean_fields(knon, knindex, tsurf_new, u0_new, v0_new)
418!
419! This routine returns the field for the ocean that has been read from the coupler
420! (done earlier with cpl_receive_frac). The field is the temperature.
421! The temperature is transformed into 1D array with valid points from index 1 to knon.
422!
423    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_ocn_day
424    USE indice_sol_mod
425
426! Input arguments
427!*************************************************************************************
428    INTEGER, INTENT(IN)                     :: knon
429    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
430
431! Output arguments
432!*************************************************************************************
433    REAL, DIMENSION(klon), INTENT(OUT)      :: tsurf_new
434    REAL, DIMENSION(klon), INTENT(OUT)      :: u0_new
435    REAL, DIMENSION(klon), INTENT(OUT)      :: v0_new
436
437! Local variables
438!*************************************************************************************
439    INTEGER                  :: i
440    INTEGER, DIMENSION(klon) :: index
441    REAL, DIMENSION(klon)    :: sic_new
442
443!*************************************************************************************
444! Transform read_sst into compressed 1D variable tsurf_new
445!
446!*************************************************************************************
447    CALL cpl2gath(read_sst, tsurf_new, knon, knindex)
448    CALL cpl2gath(read_sic, sic_new, knon, knindex)
449    CALL cpl2gath(read_u0, u0_new, knon, knindex)
450    CALL cpl2gath(read_v0, v0_new, knon, knindex)
451
452!*************************************************************************************
453! Transform read_co2 into uncompressed 1D variable fco2_ocn_day added directly in
454! the module carbon_cycle_mod
455!
456!*************************************************************************************
457    IF (carbon_cycle_cpl) THEN
458       DO i=1,klon
459          index(i)=i
460       END DO
461       CALL cpl2gath(read_co2, fco2_ocn_day, klon, index)
462    END IF
463
464!*************************************************************************************
465! The fields received from the coupler have to be weighted with the fraction of ocean
466! in relation to the total sea-ice+ocean
467!
468!*************************************************************************************
469    DO i=1, knon
470       tsurf_new(i) = tsurf_new(i)/(1. - sic_new(i))
471    END DO
472
473  END SUBROUTINE cpl_receive_ocean_fields
474
475!
476!*************************************************************************************
477!
478
479  SUBROUTINE cpl_receive_seaice_fields(knon, knindex, &
480       tsurf_new, alb_new, u0_new, v0_new)
481!
482! This routine returns the fields for the seaice that have been read from the coupler
483! (done earlier with cpl_receive_frac). These fields are the temperature and
484! albedo at sea ice surface and fraction of sea ice.
485! The fields are transformed into 1D arrays with valid points from index 1 to knon.
486!
487
488! Input arguments
489!*************************************************************************************
490    INTEGER, INTENT(IN)                     :: knon
491    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
492
493! Output arguments
494!*************************************************************************************
495    REAL, DIMENSION(klon), INTENT(OUT)      :: tsurf_new
496    REAL, DIMENSION(klon), INTENT(OUT)      :: alb_new
497    REAL, DIMENSION(klon), INTENT(OUT)      :: u0_new
498    REAL, DIMENSION(klon), INTENT(OUT)      :: v0_new
499
500! Local variables
501!*************************************************************************************
502    INTEGER               :: i
503    REAL, DIMENSION(klon) :: sic_new
504
505!*************************************************************************************
506! Transform fields read from coupler from 2D into compressed 1D variables
507!
508!*************************************************************************************
509    CALL cpl2gath(read_sit, tsurf_new, knon, knindex)
510    CALL cpl2gath(read_alb_sic, alb_new, knon, knindex)
511    CALL cpl2gath(read_sic, sic_new, knon, knindex)
512    CALL cpl2gath(read_u0, u0_new, knon, knindex)
513    CALL cpl2gath(read_v0, v0_new, knon, knindex)
514
515!*************************************************************************************
516! The fields received from the coupler have to be weighted with the sea-ice
517! concentration (in relation to the total sea-ice + ocean).
518!
519!*************************************************************************************
520    DO i= 1, knon
521       tsurf_new(i) = tsurf_new(i) / sic_new(i)
522       alb_new(i)   = alb_new(i)   / sic_new(i)
523    END DO
524
525  END SUBROUTINE cpl_receive_seaice_fields
526
527!
528!*************************************************************************************
529!
530
531  SUBROUTINE cpl_send_ocean_fields(itime, knon, knindex, &
532       swdown, lwdown, fluxlat, fluxsens, &
533       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy, windsp)
534!
535! This subroutine cumulates some fields for each time-step during a coupling
536! period. At last time-step in a coupling period the fields are transformed to the
537! grid accepted by the coupler. No sending to the coupler will be done from here
538! (it is done in cpl_send_seaice_fields).
539!
540    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
541    USE indice_sol_mod
542    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
543
544! Input arguments
545!*************************************************************************************
546    INTEGER, INTENT(IN)                     :: itime
547    INTEGER, INTENT(IN)                     :: knon
548    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
549    REAL, DIMENSION(klon), INTENT(IN)       :: swdown, lwdown
550    REAL, DIMENSION(klon), INTENT(IN)       :: fluxlat, fluxsens
551    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
552    REAL, DIMENSION(klon), INTENT(IN)       :: evap, tsurf, fder, albsol
553    REAL, DIMENSION(klon), INTENT(IN)       :: taux, tauy, windsp
554
555! Local variables
556!*************************************************************************************
557    INTEGER                                 :: cpl_index, ig
558    INTEGER                                 :: error, sum_error
559    CHARACTER(len = 25)                     :: modname = 'cpl_send_ocean_fields'
560    CHARACTER(len = 80)                     :: abort_message
561
562!*************************************************************************************
563! Start calculation
564! The ocean points are saved with second array index=1
565!
566!*************************************************************************************
567    cpl_index = 1
568
569!*************************************************************************************
570! Reset fields to zero in the beginning of a new coupling period
571!
572!*************************************************************************************
573    IF (MOD(itime, nexca) == 1) THEN
574       cpl_sols(1:knon,cpl_index) = 0.0
575       cpl_nsol(1:knon,cpl_index) = 0.0
576       cpl_rain(1:knon,cpl_index) = 0.0
577       cpl_snow(1:knon,cpl_index) = 0.0
578       cpl_evap(1:knon,cpl_index) = 0.0
579       cpl_tsol(1:knon,cpl_index) = 0.0
580       cpl_fder(1:knon,cpl_index) = 0.0
581       cpl_albe(1:knon,cpl_index) = 0.0
582       cpl_taux(1:knon,cpl_index) = 0.0
583       cpl_tauy(1:knon,cpl_index) = 0.0
584       cpl_windsp(1:knon,cpl_index) = 0.0
585       cpl_taumod(1:knon,cpl_index) = 0.0
586       IF (carbon_cycle_cpl) cpl_atm_co2(1:knon,cpl_index) = 0.0
587    ENDIF
588       
589!*************************************************************************************
590! Cumulate at each time-step
591!
592!*************************************************************************************   
593    DO ig = 1, knon
594       cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) + &
595            swdown(ig)      / REAL(nexca)
596       cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) + &
597            (lwdown(ig) + fluxlat(ig) +fluxsens(ig)) / REAL(nexca)
598       cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) + &
599            precip_rain(ig) / REAL(nexca)
600       cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) + &
601            precip_snow(ig) / REAL(nexca)
602       cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) + &
603            evap(ig)        / REAL(nexca)
604       cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) + &
605            tsurf(ig)       / REAL(nexca)
606       cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) + &
607            fder(ig)        / REAL(nexca)
608       cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) + &
609            albsol(ig)      / REAL(nexca)
610       cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) + &
611            taux(ig)        / REAL(nexca)
612       cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) + &
613            tauy(ig)        / REAL(nexca)     
614       cpl_windsp(ig,cpl_index) = cpl_windsp(ig,cpl_index) + &
615            windsp(ig)      / REAL(nexca)
616       cpl_taumod(ig,cpl_index) =   cpl_taumod(ig,cpl_index) + &
617          SQRT ( taux(ig)*taux(ig)+tauy(ig)*tauy(ig) ) / REAL (nexca)
618
619       IF (carbon_cycle_cpl) THEN
620          cpl_atm_co2(ig,cpl_index) = cpl_atm_co2(ig,cpl_index) + &
621               co2_send(knindex(ig))/ REAL(nexca)
622       END IF
623     ENDDO
624
625!*************************************************************************************
626! If the time-step corresponds to the end of coupling period the
627! fields are transformed to the 2D grid.
628! No sending to the coupler (it is done from cpl_send_seaice_fields).
629!
630!*************************************************************************************
631    IF (MOD(itime, nexca) == 0) THEN
632
633       IF (.NOT. ALLOCATED(cpl_sols2D)) THEN
634          sum_error = 0
635          ALLOCATE(cpl_sols2D(nbp_lon,jj_nb,2), stat=error)
636          sum_error = sum_error + error
637          ALLOCATE(cpl_nsol2D(nbp_lon,jj_nb,2), stat=error)
638          sum_error = sum_error + error
639          ALLOCATE(cpl_rain2D(nbp_lon,jj_nb,2), stat=error)
640          sum_error = sum_error + error
641          ALLOCATE(cpl_snow2D(nbp_lon,jj_nb,2), stat=error)
642          sum_error = sum_error + error
643          ALLOCATE(cpl_evap2D(nbp_lon,jj_nb,2), stat=error)
644          sum_error = sum_error + error
645          ALLOCATE(cpl_tsol2D(nbp_lon,jj_nb,2), stat=error)
646          sum_error = sum_error + error
647          ALLOCATE(cpl_fder2D(nbp_lon,jj_nb,2), stat=error)
648          sum_error = sum_error + error
649          ALLOCATE(cpl_albe2D(nbp_lon,jj_nb,2), stat=error)
650          sum_error = sum_error + error
651          ALLOCATE(cpl_taux2D(nbp_lon,jj_nb,2), stat=error)
652          sum_error = sum_error + error
653          ALLOCATE(cpl_tauy2D(nbp_lon,jj_nb,2), stat=error)
654          sum_error = sum_error + error
655          ALLOCATE(cpl_windsp2D(nbp_lon,jj_nb), stat=error)
656          sum_error = sum_error + error
657          ALLOCATE(cpl_taumod2D(nbp_lon,jj_nb,2), stat=error)
658          sum_error = sum_error + error
659         
660          IF (carbon_cycle_cpl) THEN
661             ALLOCATE(cpl_atm_co22D(nbp_lon,jj_nb), stat=error)
662             sum_error = sum_error + error
663          END IF
664
665          IF (sum_error /= 0) THEN
666             abort_message='Pb allocation variables couplees pour l''ecriture'
667             CALL abort_physic(modname,abort_message,1)
668          ENDIF
669       ENDIF
670       
671
672       CALL gath2cpl(cpl_sols(:,cpl_index), cpl_sols2D(:,:,cpl_index), &
673            knon, knindex)
674
675       CALL gath2cpl(cpl_nsol(:,cpl_index), cpl_nsol2D(:,:,cpl_index), &
676            knon, knindex)
677
678       CALL gath2cpl(cpl_rain(:,cpl_index), cpl_rain2D(:,:,cpl_index), &
679            knon, knindex)
680
681       CALL gath2cpl(cpl_snow(:,cpl_index), cpl_snow2D(:,:,cpl_index), &
682            knon, knindex)
683
684       CALL gath2cpl(cpl_evap(:,cpl_index), cpl_evap2D(:,:,cpl_index), &
685            knon, knindex)
686
687! cpl_tsol2D(:,:,:) not used!
688       CALL gath2cpl(cpl_tsol(:,cpl_index), cpl_tsol2D(:,:, cpl_index), &
689            knon, knindex)
690
691! cpl_fder2D(:,:,1) not used, only cpl_fder(:,:,2)!
692       CALL gath2cpl(cpl_fder(:,cpl_index), cpl_fder2D(:,:,cpl_index), &
693            knon, knindex)
694
695! cpl_albe2D(:,:,:) not used!
696       CALL gath2cpl(cpl_albe(:,cpl_index), cpl_albe2D(:,:,cpl_index), &
697            knon, knindex)
698
699       CALL gath2cpl(cpl_taux(:,cpl_index), cpl_taux2D(:,:,cpl_index), &
700            knon, knindex)
701
702       CALL gath2cpl(cpl_tauy(:,cpl_index), cpl_tauy2D(:,:,cpl_index), &
703            knon, knindex)
704
705       CALL gath2cpl(cpl_windsp(:,cpl_index), cpl_windsp2D(:,:), &
706            knon, knindex)
707
708       CALL gath2cpl(cpl_taumod(:,cpl_index), cpl_taumod2D(:,:,cpl_index), &
709            knon, knindex)
710
711       IF (carbon_cycle_cpl) &
712            CALL gath2cpl(cpl_atm_co2(:,cpl_index), cpl_atm_co22D(:,:), knon, knindex)
713   ENDIF
714
715  END SUBROUTINE cpl_send_ocean_fields
716
717!
718!*************************************************************************************
719!
720
721  SUBROUTINE cpl_send_seaice_fields(itime, dtime, knon, knindex, &
722       pctsrf, lafin, rlon, rlat, &
723       swdown, lwdown, fluxlat, fluxsens, &
724       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy)
725!
726! This subroutine cumulates some fields for each time-step during a coupling
727! period. At last time-step in a coupling period the fields are transformed to the
728! grid accepted by the coupler. All fields for all types of surfaces are sent to
729! the coupler.
730!
731    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl
732    USE indice_sol_mod
733    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
734
735! Input arguments
736!*************************************************************************************
737    INTEGER, INTENT(IN)                     :: itime
738    INTEGER, INTENT(IN)                     :: knon
739    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
740    REAL, INTENT(IN)                        :: dtime
741    REAL, DIMENSION(klon), INTENT(IN)       :: rlon, rlat
742    REAL, DIMENSION(klon), INTENT(IN)       :: swdown, lwdown
743    REAL, DIMENSION(klon), INTENT(IN)       :: fluxlat, fluxsens
744    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
745    REAL, DIMENSION(klon), INTENT(IN)       :: evap, tsurf, fder
746    REAL, DIMENSION(klon), INTENT(IN)       :: albsol, taux, tauy
747    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
748    LOGICAL, INTENT(IN)                     :: lafin
749
750! Local variables
751!*************************************************************************************
752    INTEGER                                 :: cpl_index, ig
753    INTEGER                                 :: error, sum_error
754    CHARACTER(len = 25)                     :: modname = 'cpl_send_seaice_fields'
755    CHARACTER(len = 80)                     :: abort_message
756    REAL, DIMENSION(klon)                   :: cpl_fder_tmp
757
758!*************************************************************************************
759! Start calulation
760! The sea-ice points are saved with second array index=2
761!
762!*************************************************************************************
763    cpl_index = 2
764
765!*************************************************************************************
766! Reset fields to zero in the beginning of a new coupling period
767!
768!*************************************************************************************
769    IF (MOD(itime, nexca) == 1) THEN
770       cpl_sols(1:knon,cpl_index) = 0.0
771       cpl_nsol(1:knon,cpl_index) = 0.0
772       cpl_rain(1:knon,cpl_index) = 0.0
773       cpl_snow(1:knon,cpl_index) = 0.0
774       cpl_evap(1:knon,cpl_index) = 0.0
775       cpl_tsol(1:knon,cpl_index) = 0.0
776       cpl_fder(1:knon,cpl_index) = 0.0
777       cpl_albe(1:knon,cpl_index) = 0.0
778       cpl_taux(1:knon,cpl_index) = 0.0
779       cpl_tauy(1:knon,cpl_index) = 0.0
780       cpl_taumod(1:knon,cpl_index) = 0.0
781    ENDIF
782       
783!*************************************************************************************
784! Cumulate at each time-step
785!
786!*************************************************************************************   
787    DO ig = 1, knon
788       cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) + &
789            swdown(ig)      / REAL(nexca)
790       cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) + &
791            (lwdown(ig) + fluxlat(ig) +fluxsens(ig)) / REAL(nexca)
792       cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) + &
793            precip_rain(ig) / REAL(nexca)
794       cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) + &
795            precip_snow(ig) / REAL(nexca)
796       cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) + &
797            evap(ig)        / REAL(nexca)
798       cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) + &
799            tsurf(ig)       / REAL(nexca)
800       cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) + &
801            fder(ig)        / REAL(nexca)
802       cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) + &
803            albsol(ig)      / REAL(nexca)
804       cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) + &
805            taux(ig)        / REAL(nexca)
806       cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) + &
807            tauy(ig)        / REAL(nexca)     
808       cpl_taumod(ig,cpl_index) = cpl_taumod(ig,cpl_index) + &
809            SQRT ( taux(ig)*taux(ig)+tauy(ig)*tauy(ig) ) / REAL(nexca)
810    ENDDO
811
812!*************************************************************************************
813! If the time-step corresponds to the end of coupling period the
814! fields are transformed to the 2D grid and all fields are sent to coupler.
815!
816!*************************************************************************************
817    IF (MOD(itime, nexca) == 0) THEN
818       IF (.NOT. ALLOCATED(cpl_sols2D)) THEN
819          sum_error = 0
820          ALLOCATE(cpl_sols2D(nbp_lon,jj_nb,2), stat=error)
821          sum_error = sum_error + error
822          ALLOCATE(cpl_nsol2D(nbp_lon,jj_nb,2), stat=error)
823          sum_error = sum_error + error
824          ALLOCATE(cpl_rain2D(nbp_lon,jj_nb,2), stat=error)
825          sum_error = sum_error + error
826          ALLOCATE(cpl_snow2D(nbp_lon,jj_nb,2), stat=error)
827          sum_error = sum_error + error
828          ALLOCATE(cpl_evap2D(nbp_lon,jj_nb,2), stat=error)
829          sum_error = sum_error + error
830          ALLOCATE(cpl_tsol2D(nbp_lon,jj_nb,2), stat=error)
831          sum_error = sum_error + error
832          ALLOCATE(cpl_fder2D(nbp_lon,jj_nb,2), stat=error)
833          sum_error = sum_error + error
834          ALLOCATE(cpl_albe2D(nbp_lon,jj_nb,2), stat=error)
835          sum_error = sum_error + error
836          ALLOCATE(cpl_taux2D(nbp_lon,jj_nb,2), stat=error)
837          sum_error = sum_error + error
838          ALLOCATE(cpl_tauy2D(nbp_lon,jj_nb,2), stat=error)
839          sum_error = sum_error + error
840          ALLOCATE(cpl_windsp2D(nbp_lon,jj_nb), stat=error)
841          sum_error = sum_error + error
842          ALLOCATE(cpl_taumod2D(nbp_lon,jj_nb,2), stat=error)
843          sum_error = sum_error + error
844
845          IF (carbon_cycle_cpl) THEN
846             ALLOCATE(cpl_atm_co22D(nbp_lon,jj_nb), stat=error)
847             sum_error = sum_error + error
848          END IF
849
850          IF (sum_error /= 0) THEN
851             abort_message='Pb allocation variables couplees pour l''ecriture'
852             CALL abort_physic(modname,abort_message,1)
853          ENDIF
854       ENDIF
855
856       CALL gath2cpl(cpl_sols(:,cpl_index), cpl_sols2D(:,:,cpl_index), &
857            knon, knindex)
858
859       CALL gath2cpl(cpl_nsol(:,cpl_index), cpl_nsol2D(:,:,cpl_index), &
860            knon, knindex)
861
862       CALL gath2cpl(cpl_rain(:,cpl_index), cpl_rain2D(:,:,cpl_index), &
863            knon, knindex)
864
865       CALL gath2cpl(cpl_snow(:,cpl_index), cpl_snow2D(:,:,cpl_index), &
866            knon, knindex)
867
868       CALL gath2cpl(cpl_evap(:,cpl_index), cpl_evap2D(:,:,cpl_index), &
869            knon, knindex)
870
871! cpl_tsol2D(:,:,:) not used!
872       CALL gath2cpl(cpl_tsol(:,cpl_index), cpl_tsol2D(:,:, cpl_index), &
873            knon, knindex)
874
875       ! Set default value and decompress before gath2cpl
876       cpl_fder_tmp(:) = -20.
877       DO ig = 1, knon
878          cpl_fder_tmp(knindex(ig))=cpl_fder(ig,cpl_index)
879       END DO
880       CALL gath2cpl(cpl_fder_tmp(:), cpl_fder2D(:,:,cpl_index), &
881            klon, unity)
882
883! cpl_albe2D(:,:,:) not used!
884       CALL gath2cpl(cpl_albe(:,cpl_index), cpl_albe2D(:,:,cpl_index), &
885            knon, knindex)
886
887       CALL gath2cpl(cpl_taux(:,cpl_index), cpl_taux2D(:,:,cpl_index), &
888            knon, knindex)
889
890       CALL gath2cpl(cpl_tauy(:,cpl_index), cpl_tauy2D(:,:,cpl_index), &
891            knon, knindex)
892
893       CALL gath2cpl(cpl_taumod(:,cpl_index), cpl_taumod2D(:,:,cpl_index), &
894            knon, knindex)
895
896       ! Send all fields
897       CALL cpl_send_all(itime, dtime, pctsrf, lafin, rlon, rlat)
898    ENDIF
899
900  END SUBROUTINE cpl_send_seaice_fields
901
902!
903!*************************************************************************************
904!
905
906  SUBROUTINE cpl_send_land_fields(itime, knon, knindex, rriv_in, rcoa_in)
907!
908! This subroutine cumulates some fields for each time-step during a coupling
909! period. At last time-step in a coupling period the fields are transformed to the
910! grid accepted by the coupler. No sending to the coupler will be done from here
911! (it is done in cpl_send_seaice_fields).
912!
913    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
914
915! Input arguments
916!*************************************************************************************
917    INTEGER, INTENT(IN)                       :: itime
918    INTEGER, INTENT(IN)                       :: knon
919    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex
920    REAL, DIMENSION(klon), INTENT(IN)         :: rriv_in
921    REAL, DIMENSION(klon), INTENT(IN)         :: rcoa_in
922
923! Local variables
924!*************************************************************************************
925    REAL, DIMENSION(nbp_lon,jj_nb)             :: rriv2D
926    REAL, DIMENSION(nbp_lon,jj_nb)             :: rcoa2D
927
928!*************************************************************************************
929! Rearrange fields in 2D variables
930! First initialize to zero to avoid unvalid points causing problems
931!
932!*************************************************************************************
933!$OMP MASTER
934    rriv2D(:,:) = 0.0
935    rcoa2D(:,:) = 0.0
936!$OMP END MASTER
937    CALL gath2cpl(rriv_in, rriv2D, knon, knindex)
938    CALL gath2cpl(rcoa_in, rcoa2D, knon, knindex)
939
940!*************************************************************************************
941! Reset cumulated fields to zero in the beginning of a new coupling period
942!
943!*************************************************************************************
944    IF (MOD(itime, nexca) == 1) THEN
945!$OMP MASTER
946       cpl_rriv2D(:,:) = 0.0
947       cpl_rcoa2D(:,:) = 0.0
948!$OMP END MASTER
949    ENDIF
950
951!*************************************************************************************
952! Cumulate : Following fields should be cumulated at each time-step
953!
954!*************************************************************************************   
955!$OMP MASTER
956    cpl_rriv2D(:,:) = cpl_rriv2D(:,:) + rriv2D(:,:) / REAL(nexca)
957    cpl_rcoa2D(:,:) = cpl_rcoa2D(:,:) + rcoa2D(:,:) / REAL(nexca)
958!$OMP END MASTER
959
960  END SUBROUTINE cpl_send_land_fields
961
962!
963!*************************************************************************************
964!
965
966  SUBROUTINE cpl_send_landice_fields(itime, knon, knindex, rlic_in)
967! This subroutine cumulates the field for melting ice for each time-step
968! during a coupling period. This routine will not send to coupler. Sending
969! will be done in cpl_send_seaice_fields.
970!
971
972    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
973
974! Input varibales
975!*************************************************************************************
976    INTEGER, INTENT(IN)                       :: itime
977    INTEGER, INTENT(IN)                       :: knon
978    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex
979    REAL, DIMENSION(klon), INTENT(IN)         :: rlic_in
980
981! Local varibales
982!*************************************************************************************
983    REAL, DIMENSION(nbp_lon,jj_nb)             :: rlic2D
984
985!*************************************************************************************
986! Rearrange field in a 2D variable
987! First initialize to zero to avoid unvalid points causing problems
988!
989!*************************************************************************************
990!$OMP MASTER
991    rlic2D(:,:) = 0.0
992!$OMP END MASTER
993    CALL gath2cpl(rlic_in, rlic2D, knon, knindex)
994
995!*************************************************************************************
996! Reset field to zero in the beginning of a new coupling period
997!
998!*************************************************************************************
999    IF (MOD(itime, nexca) == 1) THEN
1000!$OMP MASTER
1001       cpl_rlic2D(:,:) = 0.0
1002!$OMP END MASTER
1003    ENDIF
1004
1005!*************************************************************************************
1006! Cumulate : Melting ice should be cumulated at each time-step
1007!
1008!*************************************************************************************   
1009!$OMP MASTER
1010    cpl_rlic2D(:,:) = cpl_rlic2D(:,:) + rlic2D(:,:) / REAL(nexca)
1011!$OMP END MASTER
1012
1013  END SUBROUTINE cpl_send_landice_fields
1014
1015!
1016!*************************************************************************************
1017!
1018
1019  SUBROUTINE cpl_send_all(itime, dtime, pctsrf, lafin, rlon, rlat)
1020! This routine will send fields for all different surfaces to the coupler.
1021! This subroutine should be executed after calculations by the last surface(sea-ice),
1022! all calculations at the different surfaces have to be done before.
1023!   
1024    USE surface_data
1025    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl
1026    USE indice_sol_mod
1027    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
1028    USE time_phylmdz_mod, ONLY: start_time, itau_phy
1029! Some includes
1030!   
1031! Input arguments
1032!*************************************************************************************
1033    INTEGER, INTENT(IN)                                  :: itime
1034    REAL, INTENT(IN)                                     :: dtime
1035    REAL, DIMENSION(klon), INTENT(IN)                    :: rlon, rlat
1036    REAL, DIMENSION(klon,nbsrf), INTENT(IN)              :: pctsrf
1037    LOGICAL, INTENT(IN)                                  :: lafin
1038   
1039! Local variables
1040!*************************************************************************************
1041    INTEGER                                              :: error, sum_error, j
1042    INTEGER                                              :: itau_w
1043    INTEGER                                              :: time_sec
1044    INTEGER, DIMENSION(nbp_lon*(nbp_lat))                      :: ndexct
1045    REAL                                                 :: Up, Down
1046    REAL, DIMENSION(nbp_lon, jj_nb)                          :: tmp_lon, tmp_lat
1047    REAL, DIMENSION(nbp_lon, jj_nb, 4)                       :: pctsrf2D
1048    REAL, DIMENSION(nbp_lon, jj_nb)                          :: deno
1049    CHARACTER(len = 20)                                  :: modname = 'cpl_send_all'
1050    CHARACTER(len = 80)                                  :: abort_message
1051   
1052! Variables with fields to coupler
1053    REAL, DIMENSION(nbp_lon, jj_nb)                          :: tmp_taux
1054    REAL, DIMENSION(nbp_lon, jj_nb)                          :: tmp_tauy
1055    REAL, DIMENSION(nbp_lon, jj_nb)                          :: tmp_calv
1056! Table with all fields to send to coupler
1057    REAL, DIMENSION(nbp_lon, jj_nb, maxsend)                 :: tab_flds
1058    REAL, DIMENSION(klon_mpi)                            :: rlon_mpi, rlat_mpi
1059
1060#ifdef CPP_MPI
1061    INCLUDE 'mpif.h'
1062    INTEGER, DIMENSION(MPI_STATUS_SIZE)                  :: status
1063#endif
1064
1065! End definitions
1066!*************************************************************************************
1067   
1068
1069
1070!*************************************************************************************
1071! All fields are stored in a table tab_flds(:,:,:)
1072! First store the fields which are already on the right format
1073!
1074!*************************************************************************************
1075!$OMP MASTER
1076    tab_flds(:,:,ids_windsp) = cpl_windsp2D(:,:)
1077    tab_flds(:,:,ids_shfice) = cpl_sols2D(:,:,2)
1078    tab_flds(:,:,ids_nsfice) = cpl_nsol2D(:,:,2)
1079    tab_flds(:,:,ids_dflxdt) = cpl_fder2D(:,:,2)
1080   
1081    IF (version_ocean=='nemo') THEN
1082       tab_flds(:,:,ids_liqrun) = cpl_rriv2D(:,:) + cpl_rcoa2D(:,:)
1083       IF (carbon_cycle_cpl) tab_flds(:,:,ids_atmco2)=cpl_atm_co22D(:,:)
1084    ELSE IF (version_ocean=='opa8') THEN
1085       tab_flds(:,:,ids_shfoce) = cpl_sols2D(:,:,1)
1086       tab_flds(:,:,ids_nsfoce) = cpl_nsol2D(:,:,1)
1087       tab_flds(:,:,ids_icevap) = cpl_evap2D(:,:,2)
1088       tab_flds(:,:,ids_ocevap) = cpl_evap2D(:,:,1)
1089       tab_flds(:,:,ids_runcoa) = cpl_rcoa2D(:,:)
1090       tab_flds(:,:,ids_rivflu) = cpl_rriv2D(:,:)
1091    END IF
1092
1093!*************************************************************************************
1094! Transform the fraction of sub-surfaces from 1D to 2D array
1095!
1096!*************************************************************************************
1097    pctsrf2D(:,:,:) = 0.
1098!$OMP END MASTER
1099    CALL gath2cpl(pctsrf(:,is_oce), pctsrf2D(:,:,is_oce), klon, unity)
1100    CALL gath2cpl(pctsrf(:,is_sic), pctsrf2D(:,:,is_sic), klon, unity)
1101    CALL gath2cpl(pctsrf(:,is_lic), pctsrf2D(:,:,is_lic), klon, unity)
1102
1103!*************************************************************************************
1104! Calculate the average calving per latitude
1105! Store calving in tab_flds(:,:,19)
1106!
1107!*************************************************************************************     
1108    IF (is_omp_root) THEN
1109
1110      DO j = 1, jj_nb
1111         tmp_calv(:,j) = DOT_PRODUCT (cpl_rlic2D(1:nbp_lon,j), &
1112              pctsrf2D(1:nbp_lon,j,is_lic)) / REAL(nbp_lon)
1113      ENDDO
1114   
1115   
1116      IF (is_parallel) THEN
1117         IF (.NOT. is_north_pole_dyn) THEN
1118#ifdef CPP_MPI
1119            CALL MPI_RECV(Up,1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,status,error)
1120            CALL MPI_SEND(tmp_calv(1,1),1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,error)
1121#endif
1122         ENDIF
1123       
1124         IF (.NOT. is_south_pole_dyn) THEN
1125#ifdef CPP_MPI
1126            CALL MPI_SEND(tmp_calv(1,jj_nb),1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,error)
1127            CALL MPI_RECV(down,1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,status,error)
1128#endif
1129         ENDIF
1130         
1131         IF (.NOT. is_north_pole_dyn .AND. ii_begin /=1) THEN
1132            Up=Up+tmp_calv(nbp_lon,1)
1133            tmp_calv(:,1)=Up
1134         ENDIF
1135         
1136         IF (.NOT. is_south_pole_dyn .AND. ii_end /= nbp_lon) THEN
1137            Down=Down+tmp_calv(1,jj_nb)
1138            tmp_calv(:,jj_nb)=Down       
1139         ENDIF
1140      ENDIF
1141     
1142      tab_flds(:,:,ids_calvin) = tmp_calv(:,:)
1143
1144!*************************************************************************************
1145! Calculate total flux for snow, rain and wind with weighted addition using the
1146! fractions of ocean and seaice.
1147!
1148!*************************************************************************************   
1149       ! fraction oce+seaice
1150       deno =  pctsrf2D(:,:,is_oce) + pctsrf2D(:,:,is_sic)
1151
1152       IF (version_ocean=='nemo') THEN
1153          tab_flds(:,:,ids_shftot)  = 0.0
1154          tab_flds(:,:,ids_nsftot) = 0.0
1155          tab_flds(:,:,ids_totrai) = 0.0
1156          tab_flds(:,:,ids_totsno) = 0.0
1157          tab_flds(:,:,ids_toteva) = 0.0
1158          tab_flds(:,:,ids_taumod) = 0.0
1159 
1160          tmp_taux(:,:)    = 0.0
1161          tmp_tauy(:,:)    = 0.0
1162          ! For all valid grid cells containing some fraction of ocean or sea-ice
1163          WHERE ( deno(:,:) /= 0 )
1164             tmp_taux = cpl_taux2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1165                  cpl_taux2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1166             tmp_tauy = cpl_tauy2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1167                  cpl_tauy2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1168
1169             tab_flds(:,:,ids_shftot) = cpl_sols2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1170                  cpl_sols2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1171             tab_flds(:,:,ids_nsftot) = cpl_nsol2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1172                  cpl_nsol2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1173             tab_flds(:,:,ids_totrai) = cpl_rain2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1174                  cpl_rain2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1175             tab_flds(:,:,ids_totsno) = cpl_snow2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1176                  cpl_snow2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1177             tab_flds(:,:,ids_toteva) = cpl_evap2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1178                  cpl_evap2D(:,:,2)  * pctsrf2D(:,:,is_sic) / deno(:,:)
1179             tab_flds(:,:,ids_taumod) = cpl_taumod2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1180                  cpl_taumod2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1181             
1182         ENDWHERE
1183
1184          tab_flds(:,:,ids_icevap) = cpl_evap2D(:,:,2)
1185         
1186       ELSE IF (version_ocean=='opa8') THEN
1187          ! Store fields for rain and snow in tab_flds(:,:,15) and tab_flds(:,:,16)
1188          tab_flds(:,:,ids_totrai) = 0.0
1189          tab_flds(:,:,ids_totsno) = 0.0
1190          tmp_taux(:,:)    = 0.0
1191          tmp_tauy(:,:)    = 0.0
1192          ! For all valid grid cells containing some fraction of ocean or sea-ice
1193          WHERE ( deno(:,:) /= 0 )
1194             tab_flds(:,:,ids_totrai) = cpl_rain2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1195                  cpl_rain2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1196             tab_flds(:,:,ids_totsno) = cpl_snow2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1197                  cpl_snow2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1198             
1199             tmp_taux = cpl_taux2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1200                  cpl_taux2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1201             tmp_tauy = cpl_tauy2D(:,:,1) * pctsrf2D(:,:,is_oce) / deno(:,:) +    &
1202                  cpl_tauy2D(:,:,2) * pctsrf2D(:,:,is_sic) / deno(:,:)
1203          ENDWHERE
1204       END IF
1205
1206    ENDIF ! is_omp_root
1207 
1208!*************************************************************************************
1209! Transform the wind components from local atmospheric 2D coordinates to geocentric
1210! 3D coordinates.
1211! Store the resulting wind components in tab_flds(:,:,1:6)
1212!*************************************************************************************
1213
1214! Transform the longitudes and latitudes on 2D arrays
1215   
1216    CALL gather_omp(rlon,rlon_mpi)
1217    CALL gather_omp(rlat,rlat_mpi)
1218!$OMP MASTER
1219    CALL Grid1DTo2D_mpi(rlon_mpi,tmp_lon)
1220    CALL Grid1DTo2D_mpi(rlat_mpi,tmp_lat)
1221!$OMP END MASTER   
1222
1223    IF (is_sequential) THEN
1224       IF (is_north_pole_dyn) tmp_lon(:,1)     = tmp_lon(:,2)
1225       IF (is_south_pole_dyn) tmp_lon(:,nbp_lat) = tmp_lon(:,nbp_lat-1)
1226    ENDIF
1227     
1228! NetCDF output of the wind before transformation of coordinate system
1229    IF (is_sequential) THEN
1230       ndexct(:) = 0
1231       itau_w = itau_phy + itime + start_time * day_step_phy
1232       CALL histwrite(nidct,'tauxe',itau_w,tmp_taux,nbp_lon*(nbp_lat),ndexct)
1233       CALL histwrite(nidct,'tauyn',itau_w,tmp_tauy,nbp_lon*(nbp_lat),ndexct)
1234       CALL histwrite(nidct,'tmp_lon',itau_w,tmp_lon,nbp_lon*(nbp_lat),ndexct)
1235       CALL histwrite(nidct,'tmp_lat',itau_w,tmp_lat,nbp_lon*(nbp_lat),ndexct)
1236    ENDIF
1237
1238! Transform the wind from spherical atmospheric 2D coordinates to geocentric
1239! cartesian 3D coordinates
1240!$OMP MASTER
1241    CALL atm2geo (nbp_lon, jj_nb, tmp_taux, tmp_tauy, tmp_lon, tmp_lat, &
1242         tab_flds(:,:,ids_tauxxu), tab_flds(:,:,ids_tauyyu), tab_flds(:,:,ids_tauzzu) )
1243   
1244    tab_flds(:,:,ids_tauxxv)  = tab_flds(:,:,ids_tauxxu)
1245    tab_flds(:,:,ids_tauyyv)  = tab_flds(:,:,ids_tauyyu)
1246    tab_flds(:,:,ids_tauzzv)  = tab_flds(:,:,ids_tauzzu)
1247!$OMP END MASTER
1248
1249!*************************************************************************************
1250! NetCDF output of all fields just before sending to coupler.
1251!
1252!*************************************************************************************
1253    IF (is_sequential) THEN
1254        DO j=1,maxsend
1255          IF (infosend(j)%action) CALL histwrite(nidct,infosend(j)%name, itau_w, &
1256             tab_flds(:,:,j),nbp_lon*(nbp_lat),ndexct)
1257        ENDDO
1258    ENDIF
1259!*************************************************************************************
1260! Send the table of all fields
1261!
1262!*************************************************************************************
1263    time_sec=(itime-1)*dtime
1264#ifdef CPP_COUPLE
1265!$OMP MASTER
1266    CALL intocpl(time_sec, lafin, tab_flds(:,:,:))
1267!$OMP END MASTER
1268#endif
1269
1270!*************************************************************************************
1271! Finish with some dellocate
1272!
1273!************************************************************************************* 
1274    sum_error=0
1275    DEALLOCATE(cpl_sols2D, cpl_nsol2D, cpl_rain2D, cpl_snow2D, stat=error )
1276    sum_error = sum_error + error
1277    DEALLOCATE(cpl_evap2D, cpl_tsol2D, cpl_fder2D, cpl_albe2D, stat=error )
1278    sum_error = sum_error + error
1279    DEALLOCATE(cpl_taux2D, cpl_tauy2D, cpl_windsp2D, cpl_taumod2D, stat=error )
1280    sum_error = sum_error + error
1281   
1282    IF (carbon_cycle_cpl) THEN
1283       DEALLOCATE(cpl_atm_co22D, stat=error )
1284       sum_error = sum_error + error
1285    END IF
1286
1287    IF (sum_error /= 0) THEN
1288       abort_message='Pb in deallocation of cpl_xxxx2D coupling variables'
1289       CALL abort_physic(modname,abort_message,1)
1290    ENDIF
1291   
1292  END SUBROUTINE cpl_send_all
1293!
1294!*************************************************************************************
1295!
1296  SUBROUTINE cpl2gath(champ_in, champ_out, knon, knindex)
1297  USE mod_phys_lmdz_para
1298! Cette routine transforme un champs de la grille 2D recu du coupleur sur la grille
1299! 'gathered' (la grille physiq comprime).
1300!
1301!
1302! input:         
1303!   champ_in     champ sur la grille 2D
1304!   knon         nombre de points dans le domaine a traiter
1305!   knindex      index des points de la surface a traiter
1306!
1307! output:
1308!   champ_out    champ sur la grille 'gatherd'
1309!
1310    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
1311
1312! Input
1313    INTEGER, INTENT(IN)                       :: knon
1314    REAL, DIMENSION(nbp_lon,jj_nb), INTENT(IN)    :: champ_in
1315    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex
1316
1317! Output
1318    REAL, DIMENSION(klon_mpi), INTENT(OUT)        :: champ_out
1319
1320! Local
1321    INTEGER                                   :: i, ig
1322    REAL, DIMENSION(klon_mpi)                 :: temp_mpi
1323    REAL, DIMENSION(klon)                     :: temp_omp
1324
1325!*************************************************************************************
1326!
1327   
1328
1329! Transform from 2 dimensions (nbp_lon,jj_nb) to 1 dimension (klon)
1330!$OMP MASTER
1331    CALL Grid2Dto1D_mpi(champ_in,temp_mpi)
1332!$OMP END MASTER
1333
1334    CALL scatter_omp(temp_mpi,temp_omp)
1335   
1336! Compress from klon to knon
1337    DO i = 1, knon
1338       ig = knindex(i)
1339       champ_out(i) = temp_omp(ig)
1340    ENDDO
1341
1342  END SUBROUTINE cpl2gath
1343!
1344!*************************************************************************************
1345!
1346  SUBROUTINE gath2cpl(champ_in, champ_out, knon, knindex)
1347  USE mod_phys_lmdz_para
1348! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
1349! au coupleur.
1350!
1351! input:         
1352!   champ_in     champ sur la grille gathere       
1353!   knon         nombre de points dans le domaine a traiter
1354!   knindex      index des points de la surface a traiter
1355!
1356! output:
1357!   champ_out    champ sur la grille 2D
1358!
1359    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
1360   
1361! Input arguments
1362!*************************************************************************************
1363    INTEGER, INTENT(IN)                    :: knon
1364    REAL, DIMENSION(klon), INTENT(IN)      :: champ_in
1365    INTEGER, DIMENSION(klon), INTENT(IN)   :: knindex
1366
1367! Output arguments
1368!*************************************************************************************
1369    REAL, DIMENSION(nbp_lon,jj_nb), INTENT(OUT) :: champ_out
1370
1371! Local variables
1372!*************************************************************************************
1373    INTEGER                                :: i, ig
1374    REAL, DIMENSION(klon)                  :: temp_omp
1375    REAL, DIMENSION(klon_mpi)              :: temp_mpi
1376!*************************************************************************************
1377
1378! Decompress from knon to klon
1379    temp_omp = 0.
1380    DO i = 1, knon
1381       ig = knindex(i)
1382       temp_omp(ig) = champ_in(i)
1383    ENDDO
1384
1385! Transform from 1 dimension (klon) to 2 dimensions (nbp_lon,jj_nb)
1386    CALL gather_omp(temp_omp,temp_mpi)
1387
1388!$OMP MASTER   
1389    CALL Grid1Dto2D_mpi(temp_mpi,champ_out)
1390   
1391    IF (is_north_pole_dyn) champ_out(:,1)=temp_mpi(1)
1392    IF (is_south_pole_dyn) champ_out(:,jj_nb)=temp_mpi(klon)
1393!$OMP END MASTER
1394   
1395  END SUBROUTINE gath2cpl
1396!
1397!*************************************************************************************
1398!
1399END MODULE cpl_mod
1400
Note: See TracBrowser for help on using the repository browser.