GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/surf_land_orchidee_mod.F90 Lines: 0 284 0.0 %
Date: 2023-06-30 12:56:34 Branches: 0 416 0.0 %

Line Branch Exec Source
1
!
2
MODULE surf_land_orchidee_mod
3
#ifndef ORCHIDEE_NOOPENMP
4
#ifndef ORCHIDEE_NOZ0H
5
#ifndef ORCHIDEE_NOFREIN
6
#ifndef ORCHIDEE_NOUNSTRUCT
7
#ifndef ORCHIDEE_NOLIC
8
!
9
! This module controles the interface towards the model ORCHIDEE.
10
!
11
! Compatibility with ORCHIDIEE :
12
! The current version can be used with ORCHIDEE/trunk from revision 7757.
13
! This interface is used if none of the cpp keys ORCHIDEE_NOOPENMP,
14
! ORCHIDEE_NOZ0H, ORCHIDEE_NOFREIN or ORCHIDEE_NOLIC is set.
15
!
16
! Subroutines in this module : surf_land_orchidee
17
!                              Init_orchidee_index
18
!                              Get_orchidee_communicator
19
!                              Init_neighbours
20
21
  USE dimphy
22
#ifdef CPP_VEGET
23
  USE intersurf     ! module in ORCHIDEE
24
#endif
25
  USE cpl_mod,      ONLY : cpl_send_land_fields, cpl_send_landice_fields
26
  USE surface_data, ONLY : type_ocean, landice_opt
27
  USE geometry_mod, ONLY : dx, dy, boundslon, boundslat,longitude, latitude, cell_area,  ind_cell_glo
28
  USE mod_grid_phy_lmdz
29
  USE mod_phys_lmdz_para, mpi_root_rank=>mpi_master
30
  USE carbon_cycle_mod, ONLY : nbcf_in_orc, nbcf_out, fields_in, yfields_in, yfields_out, cfname_in, cfname_out
31
  USE nrtype, ONLY : PI
32
33
  IMPLICIT NONE
34
35
  PRIVATE
36
  PUBLIC  :: surf_land_orchidee
37
38
CONTAINS
39
!
40
!****************************************************************************************
41
!
42
  SUBROUTINE surf_land_orchidee(itime, dtime, date0, knon, &
43
       knindex, rlon, rlat, yrmu0, pctsrf, &
44
       debut, lafin, &
45
       plev,  u1_lay, v1_lay, gustiness, temp_air, spechum, epot_air, ccanopy, &
46
       tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
47
       precip_rain, precip_snow, lwdown, swnet, swdown, &
48
       ps, q2m, t2m, &
49
       evap, fluxsens, fluxlat, &
50
       tsol_rad, tsurf_new, alb1_new, alb2_new, &
51
       emis_new, z0m_new, z0h_new, qsurf, &
52
       veget, lai, height )
53
54
    USE mod_surf_para
55
    USE mod_synchro_omp
56
    USE carbon_cycle_mod
57
    USE indice_sol_mod
58
    USE print_control_mod, ONLY: lunout
59
    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
60
#ifdef CPP_VEGET
61
    USE time_phylmdz_mod, ONLY: itau_phy
62
#endif
63
!
64
! Cette routine sert d'interface entre le modele atmospherique et le
65
! modele de sol continental. Appel a sechiba
66
!
67
! L. Fairhead 02/2000
68
!
69
! input:
70
!   itime        numero du pas de temps
71
!   dtime        pas de temps de la physique (en s)
72
!   nisurf       index de la surface a traiter (1 = sol continental)
73
!   knon         nombre de points de la surface a traiter
74
!   knindex      index des points de la surface a traiter
75
!   rlon         longitudes de la grille entiere
76
!   rlat         latitudes de la grille entiere
77
!   pctsrf       tableau des fractions de surface de chaque maille
78
!   debut        logical: 1er appel a la physique (lire les restart)
79
!   lafin        logical: dernier appel a la physique (ecrire les restart)
80
!                     (si false calcul simplifie des fluxs sur les continents)
81
!   plev         hauteur de la premiere couche (Pa)
82
!   u1_lay       vitesse u 1ere couche
83
!   v1_lay       vitesse v 1ere couche
84
!   temp_air     temperature de l'air 1ere couche
85
!   spechum      humidite specifique 1ere couche
86
!   epot_air     temp pot de l'air
87
!   ccanopy      concentration CO2 canopee, correspond au co2_send de
88
!                carbon_cycle_mod ou valeur constant co2_ppm
89
!   tq_cdrag     cdrag
90
!   petAcoef     coeff. A de la resolution de la CL pour t
91
!   peqAcoef     coeff. A de la resolution de la CL pour q
92
!   petBcoef     coeff. B de la resolution de la CL pour t
93
!   peqBcoef     coeff. B de la resolution de la CL pour q
94
!   precip_rain  precipitation liquide
95
!   precip_snow  precipitation solide
96
!   lwdown       flux IR descendant a la surface
97
!   swnet        flux solaire net
98
!   swdown       flux solaire entrant a la surface
99
!   ps           pression au sol
100
!   radsol       rayonnement net aus sol (LW + SW)
101
!
102
! output:
103
!   evap         evaporation totale
104
!   fluxsens     flux de chaleur sensible
105
!   fluxlat      flux de chaleur latente
106
!   tsol_rad
107
!   tsurf_new    temperature au sol
108
!   alb1_new     albedo in visible SW interval
109
!   alb2_new     albedo in near IR interval
110
!   emis_new     emissivite
111
!   z0m_new      surface roughness for momentum
112
!   z0h_new      surface roughness for heat
113
!   qsurf        air moisture at surface
114
!
115
    INCLUDE "YOMCST.h"
116
    INCLUDE "dimpft.h"
117
!
118
! Parametres d'entree
119
!****************************************************************************************
120
    INTEGER, INTENT(IN)                       :: itime
121
    REAL, INTENT(IN)                          :: dtime
122
    REAL, INTENT(IN)                          :: date0
123
    INTEGER, INTENT(IN)                       :: knon
124
    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex
125
    LOGICAL, INTENT(IN)                       :: debut, lafin
126
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)   :: pctsrf
127
    REAL, DIMENSION(klon), INTENT(IN)         :: rlon, rlat
128
    REAL, DIMENSION(klon), INTENT(IN)         :: yrmu0 ! cosine of solar zenith angle
129
    REAL, DIMENSION(klon), INTENT(IN)         :: plev
130
    REAL, DIMENSION(klon), INTENT(IN)         :: u1_lay, v1_lay, gustiness
131
    REAL, DIMENSION(klon), INTENT(IN)         :: temp_air, spechum
132
    REAL, DIMENSION(klon), INTENT(IN)         :: epot_air, ccanopy
133
    REAL, DIMENSION(klon), INTENT(IN)         :: tq_cdrag
134
    REAL, DIMENSION(klon), INTENT(IN)         :: petAcoef, peqAcoef
135
    REAL, DIMENSION(klon), INTENT(IN)         :: petBcoef, peqBcoef
136
    REAL, DIMENSION(klon), INTENT(IN)         :: precip_rain, precip_snow
137
    REAL, DIMENSION(klon), INTENT(IN)         :: lwdown, swnet, swdown, ps
138
    REAL, DIMENSION(klon), INTENT(IN)         :: q2m, t2m
139
140
! Parametres de sortie
141
!****************************************************************************************
142
    REAL, DIMENSION(klon), INTENT(OUT)        :: evap, fluxsens, fluxlat, qsurf
143
    REAL, DIMENSION(klon), INTENT(OUT)        :: tsol_rad, tsurf_new
144
    REAL, DIMENSION(klon), INTENT(OUT)        :: alb1_new, alb2_new
145
    REAL, DIMENSION(klon), INTENT(OUT)        :: emis_new, z0m_new, z0h_new
146
    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: veget
147
    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: lai
148
    REAL, DIMENSION(klon,nvm_lmdz), INTENT(OUT) :: height
149
150
! Local
151
!****************************************************************************************
152
    INTEGER                                   :: ij, jj, igrid, ireal, index, nb
153
    INTEGER                                   :: error
154
    REAL, DIMENSION(klon)                     :: swdown_vrai
155
    REAL, DIMENSION(klon)                     :: run_off_lic        !! run off from land ice defined in ORCHIDEE, contains calving, melting and liquid precipitation
156
    REAL, DIMENSION(klon)                     :: run_off_lic_frac   !! cell fraction corresponding to run_off_lic
157
    REAL, DIMENSION(klon)                     :: blowingsnow_flux   !! blowing snow flux
158
    CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
159
    CHARACTER (len = 80)                      :: abort_message
160
    LOGICAL,SAVE                              :: check = .FALSE.
161
    !$OMP THREADPRIVATE(check)
162
163
! type de couplage dans sechiba
164
!  character (len=10)   :: coupling = 'implicit'
165
! drapeaux controlant les appels dans SECHIBA
166
!  type(control_type), save   :: control_in
167
! Preserved albedo
168
    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: albedo_keep, zlev
169
    !$OMP THREADPRIVATE(albedo_keep,zlev)
170
! coordonnees geographiques
171
    REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: lalo
172
    !$OMP THREADPRIVATE(lalo)
173
! boundaries of cells
174
    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: bounds_lalo
175
    !$OMP THREADPRIVATE(bounds_lalo)
176
! pts voisins
177
    INTEGER,ALLOCATABLE, DIMENSION(:,:), SAVE :: neighbours
178
    !$OMP THREADPRIVATE(neighbours)
179
! fractions continents
180
    REAL,ALLOCATABLE, DIMENSION(:), SAVE      :: contfrac
181
    !$OMP THREADPRIVATE(contfrac)
182
! resolution de la grille
183
    REAL, ALLOCATABLE, DIMENSION (:,:), SAVE  :: resolution
184
    !$OMP THREADPRIVATE(resolution)
185
186
    REAL, ALLOCATABLE, DIMENSION (:,:), SAVE  :: lon_scat, lat_scat
187
    !$OMP THREADPRIVATE(lon_scat,lat_scat)
188
189
! area of cells
190
    REAL, ALLOCATABLE, DIMENSION (:), SAVE  :: area
191
    !$OMP THREADPRIVATE(area)
192
193
    LOGICAL, SAVE                             :: lrestart_read = .TRUE.
194
    !$OMP THREADPRIVATE(lrestart_read)
195
    LOGICAL, SAVE                             :: lrestart_write = .FALSE.
196
    !$OMP THREADPRIVATE(lrestart_write)
197
198
    REAL, DIMENSION(knon,2)                   :: albedo_out
199
200
! Pb de nomenclature
201
    REAL, DIMENSION(klon)                     :: petA_orc, peqA_orc
202
    REAL, DIMENSION(klon)                     :: petB_orc, peqB_orc
203
! Pb de correspondances de grilles
204
    INTEGER, DIMENSION(:), SAVE, ALLOCATABLE  :: ig, jg
205
    !$OMP THREADPRIVATE(ig,jg)
206
    INTEGER :: indi, indj
207
    INTEGER, SAVE, ALLOCATABLE,DIMENSION(:)   :: ktindex
208
    !$OMP THREADPRIVATE(ktindex)
209
210
! Essai cdrag
211
    REAL, DIMENSION(klon)                     :: cdrag
212
    INTEGER,SAVE                              :: offset
213
    !$OMP THREADPRIVATE(offset)
214
215
    REAL, DIMENSION(klon_glo)                 :: rlon_g,rlat_g
216
    INTEGER, SAVE                             :: orch_comm
217
    !$OMP THREADPRIVATE(orch_comm)
218
219
    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: coastalflow
220
    !$OMP THREADPRIVATE(coastalflow)
221
    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: riverflow
222
    !$OMP THREADPRIVATE(riverflow)
223
224
    INTEGER :: orch_mpi_rank
225
    INTEGER :: orch_mpi_size
226
    INTEGER :: orch_omp_rank
227
    INTEGER :: orch_omp_size
228
229
    REAL, ALLOCATABLE, DIMENSION(:)         :: longitude_glo
230
    REAL, ALLOCATABLE, DIMENSION(:)         :: latitude_glo
231
    REAL, ALLOCATABLE, DIMENSION(:,:)       :: boundslon_glo
232
    REAL, ALLOCATABLE, DIMENSION(:,:)       :: boundslat_glo
233
    INTEGER, ALLOCATABLE, DIMENSION(:)      :: ind_cell_glo_glo
234
    INTEGER, ALLOCATABLE, SAVE,DIMENSION(:) :: ind_cell
235
    !$OMP THREADPRIVATE(ind_cell)
236
    INTEGER :: begin, end
237
!
238
! Fin definition
239
!****************************************************************************************
240
241
    IF (check) WRITE(lunout,*)'Entree ', modname
242
243
! Initialisation
244
245
    IF (debut) THEN
246
! Test of coherence between variable ok_veget and cpp key CPP_VEGET
247
#ifndef CPP_VEGET
248
       abort_message='Pb de coherence: ok_veget = .true. mais CPP_VEGET = .false.'
249
       CALL abort_physic(modname,abort_message,1)
250
#endif
251
252
       CALL Init_surf_para(knon)
253
       ALLOCATE(ktindex(knon))
254
       IF ( .NOT. ALLOCATED(albedo_keep)) THEN
255
!ym          ALLOCATE(albedo_keep(klon))
256
!ym bizarre que non allou� en knon precedement
257
          ALLOCATE(albedo_keep(knon))
258
          ALLOCATE(zlev(knon))
259
       ENDIF
260
! Pb de correspondances de grilles
261
       ALLOCATE(ig(klon))
262
       ALLOCATE(jg(klon))
263
       ig(1) = 1
264
       jg(1) = 1
265
       indi = 0
266
       indj = 2
267
       DO igrid = 2, klon - 1
268
          indi = indi + 1
269
          IF ( indi > nbp_lon) THEN
270
             indi = 1
271
             indj = indj + 1
272
          ENDIF
273
          ig(igrid) = indi
274
          jg(igrid) = indj
275
       ENDDO
276
       ig(klon) = 1
277
       jg(klon) = nbp_lat
278
279
       IF ((.NOT. ALLOCATED(area))) THEN
280
          ALLOCATE(area(knon), stat = error)
281
          IF (error /= 0) THEN
282
             abort_message='Pb allocation area'
283
             CALL abort_physic(modname,abort_message,1)
284
          ENDIF
285
       ENDIF
286
       DO igrid = 1, knon
287
          area(igrid) = cell_area(knindex(igrid))
288
       ENDDO
289
290
       IF (grid_type==unstructured) THEN
291
292
293
         IF ((.NOT. ALLOCATED(lon_scat))) THEN
294
            ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
295
            IF (error /= 0) THEN
296
               abort_message='Pb allocation lon_scat'
297
               CALL abort_physic(modname,abort_message,1)
298
            ENDIF
299
         ENDIF
300
301
         IF ((.NOT. ALLOCATED(lat_scat))) THEN
302
            ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
303
            IF (error /= 0) THEN
304
               abort_message='Pb allocation lat_scat'
305
               CALL abort_physic(modname,abort_message,1)
306
            ENDIF
307
         ENDIF
308
         CALL Gather(rlon,rlon_g)
309
         CALL Gather(rlat,rlat_g)
310
311
         IF (is_mpi_root) THEN
312
            index = 1
313
            DO jj = 2, nbp_lat-1
314
               DO ij = 1, nbp_lon
315
                  index = index + 1
316
                  lon_scat(ij,jj) = rlon_g(index)
317
                  lat_scat(ij,jj) = rlat_g(index)
318
               ENDDO
319
            ENDDO
320
            lon_scat(:,1) = lon_scat(:,2)
321
            lat_scat(:,1) = rlat_g(1)
322
            lon_scat(:,nbp_lat) = lon_scat(:,2)
323
            lat_scat(:,nbp_lat) = rlat_g(klon_glo)
324
         ENDIF
325
326
         CALL bcast(lon_scat)
327
         CALL bcast(lat_scat)
328
329
       ELSE IF (grid_type==regular_lonlat) THEN
330
331
         IF ((.NOT. ALLOCATED(lalo))) THEN
332
            ALLOCATE(lalo(knon,2), stat = error)
333
            IF (error /= 0) THEN
334
               abort_message='Pb allocation lalo'
335
               CALL abort_physic(modname,abort_message,1)
336
            ENDIF
337
         ENDIF
338
339
         IF ((.NOT. ALLOCATED(bounds_lalo))) THEN
340
           ALLOCATE(bounds_lalo(knon,nvertex,2), stat = error)
341
           IF (error /= 0) THEN
342
             abort_message='Pb allocation lalo'
343
             CALL abort_physic(modname,abort_message,1)
344
           ENDIF
345
         ENDIF
346
347
         IF ((.NOT. ALLOCATED(lon_scat))) THEN
348
            ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
349
            IF (error /= 0) THEN
350
               abort_message='Pb allocation lon_scat'
351
               CALL abort_physic(modname,abort_message,1)
352
            ENDIF
353
         ENDIF
354
         IF ((.NOT. ALLOCATED(lat_scat))) THEN
355
            ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
356
            IF (error /= 0) THEN
357
               abort_message='Pb allocation lat_scat'
358
               CALL abort_physic(modname,abort_message,1)
359
            ENDIF
360
         ENDIF
361
         lon_scat = 0.
362
         lat_scat = 0.
363
         DO igrid = 1, knon
364
            index = knindex(igrid)
365
            lalo(igrid,2) = rlon(index)
366
            lalo(igrid,1) = rlat(index)
367
            bounds_lalo(igrid,:,2)=boundslon(index,:)*180./PI
368
            bounds_lalo(igrid,:,1)=boundslat(index,:)*180./PI
369
         ENDDO
370
371
372
373
         CALL Gather(rlon,rlon_g)
374
         CALL Gather(rlat,rlat_g)
375
376
         IF (is_mpi_root) THEN
377
            index = 1
378
            DO jj = 2, nbp_lat-1
379
               DO ij = 1, nbp_lon
380
                  index = index + 1
381
                  lon_scat(ij,jj) = rlon_g(index)
382
                  lat_scat(ij,jj) = rlat_g(index)
383
               ENDDO
384
            ENDDO
385
            lon_scat(:,1) = lon_scat(:,2)
386
            lat_scat(:,1) = rlat_g(1)
387
            lon_scat(:,nbp_lat) = lon_scat(:,2)
388
            lat_scat(:,nbp_lat) = rlat_g(klon_glo)
389
         ENDIF
390
391
         CALL bcast(lon_scat)
392
         CALL bcast(lat_scat)
393
394
       ENDIF
395
!
396
! Allouer et initialiser le tableau des voisins et des fraction de continents
397
!
398
       IF (( .NOT. ALLOCATED(contfrac))) THEN
399
          ALLOCATE(contfrac(knon), stat = error)
400
          IF (error /= 0) THEN
401
             abort_message='Pb allocation contfrac'
402
             CALL abort_physic(modname,abort_message,1)
403
          ENDIF
404
       ENDIF
405
406
       DO igrid = 1, knon
407
          ireal = knindex(igrid)
408
          contfrac(igrid) = pctsrf(ireal,is_ter)
409
       ENDDO
410
411
412
       IF (grid_type==regular_lonlat) THEN
413
414
         IF ( (.NOT.ALLOCATED(neighbours))) THEN
415
          ALLOCATE(neighbours(knon,8), stat = error)
416
          IF (error /= 0) THEN
417
             abort_message='Pb allocation neighbours'
418
             CALL abort_physic(modname,abort_message,1)
419
          ENDIF
420
         ENDIF
421
         neighbours = -1.
422
         CALL Init_neighbours(knon,neighbours,knindex,pctsrf(:,is_ter))
423
424
       ELSE IF (grid_type==unstructured) THEN
425
426
         IF ( (.NOT.ALLOCATED(neighbours))) THEN
427
          ALLOCATE(neighbours(knon,12), stat = error)
428
          IF (error /= 0) THEN
429
             abort_message='Pb allocation neighbours'
430
             CALL abort_physic(modname,abort_message,1)
431
          ENDIF
432
         ENDIF
433
         neighbours = -1.
434
435
       ENDIF
436
437
438
!
439
!  Allocation et calcul resolutions
440
       IF ( (.NOT.ALLOCATED(resolution))) THEN
441
          ALLOCATE(resolution(knon,2), stat = error)
442
          IF (error /= 0) THEN
443
             abort_message='Pb allocation resolution'
444
             CALL abort_physic(modname,abort_message,1)
445
          ENDIF
446
       ENDIF
447
448
       IF (grid_type==regular_lonlat) THEN
449
         DO igrid = 1, knon
450
            ij = knindex(igrid)
451
            resolution(igrid,1) = dx(ij)
452
           resolution(igrid,2) = dy(ij)
453
         ENDDO
454
       ENDIF
455
456
       ALLOCATE(coastalflow(klon), stat = error)
457
       IF (error /= 0) THEN
458
          abort_message='Pb allocation coastalflow'
459
          CALL abort_physic(modname,abort_message,1)
460
       ENDIF
461
462
       ALLOCATE(riverflow(klon), stat = error)
463
       IF (error /= 0) THEN
464
          abort_message='Pb allocation riverflow'
465
          CALL abort_physic(modname,abort_message,1)
466
       ENDIF
467
!
468
! carbon_cycle_cpl not possible with this interface and version of ORHCHIDEE
469
!
470
! >> PC
471
!       IF (carbon_cycle_cpl) THEN
472
!          abort_message='carbon_cycle_cpl not yet possible with this interface of ORCHIDEE'
473
!          CALL abort_physic(modname,abort_message,1)
474
!       END IF
475
! << PC
476
477
    ENDIF                          ! (fin debut)
478
479
!
480
! Appel a la routine sols continentaux
481
!
482
    IF (lafin) lrestart_write = .TRUE.
483
    IF (check) WRITE(lunout,*)'lafin ',lafin,lrestart_write
484
485
    petA_orc(1:knon) = petBcoef(1:knon) * dtime
486
    petB_orc(1:knon) = petAcoef(1:knon)
487
    peqA_orc(1:knon) = peqBcoef(1:knon) * dtime
488
    peqB_orc(1:knon) = peqAcoef(1:knon)
489
490
    cdrag = 0.
491
    cdrag(1:knon) = tq_cdrag(1:knon)
492
493
! zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/287.05*temp_air(1:knon))*9.80665)
494
!    zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/RD*temp_air(1:knon))*RG)
495
     zlev(1:knon) = plev(1:knon)*RD*temp_air(1:knon)/((ps(1:knon)*100.0)*RG)
496
497
498
! PF et PASB
499
!   where(cdrag > 0.01)
500
!     cdrag = 0.01
501
!   endwhere
502
!  write(*,*)'Cdrag = ',minval(cdrag),maxval(cdrag)
503
504
505
    IF (debut) THEN
506
       CALL Init_orchidee_index(knon,knindex,offset,ktindex)
507
       CALL Get_orchidee_communicator(orch_comm,orch_mpi_size,orch_mpi_rank, orch_omp_size,orch_omp_rank)
508
509
       IF (grid_type==unstructured) THEN
510
         IF (knon==0) THEN
511
           begin=1
512
           end=0
513
         ELSE
514
           begin=offset+1
515
           end=offset+ktindex(knon)
516
         ENDIF
517
518
         IF (orch_mpi_rank==orch_mpi_size-1 .AND. orch_omp_rank==orch_omp_size-1) end=nbp_lon*nbp_lat
519
520
         ALLOCATE(lalo(end-begin+1,2))
521
         ALLOCATE(bounds_lalo(end-begin+1,nvertex,2))
522
         ALLOCATE(ind_cell(end-begin+1))
523
524
         ALLOCATE(longitude_glo(klon_glo))
525
         CALL gather(longitude,longitude_glo)
526
         CALL bcast(longitude_glo)
527
         lalo(:,2)=longitude_glo(begin:end)*180./PI
528
529
         ALLOCATE(latitude_glo(klon_glo))
530
         CALL gather(latitude,latitude_glo)
531
         CALL bcast(latitude_glo)
532
         lalo(:,1)=latitude_glo(begin:end)*180./PI
533
534
         ALLOCATE(boundslon_glo(klon_glo,nvertex))
535
         CALL gather(boundslon,boundslon_glo)
536
         CALL bcast(boundslon_glo)
537
         bounds_lalo(:,:,2)=boundslon_glo(begin:end,:)*180./PI
538
539
         ALLOCATE(boundslat_glo(klon_glo,nvertex))
540
         CALL gather(boundslat,boundslat_glo)
541
         CALL bcast(boundslat_glo)
542
         bounds_lalo(:,:,1)=boundslat_glo(begin:end,:)*180./PI
543
544
         ALLOCATE(ind_cell_glo_glo(klon_glo))
545
         CALL gather(ind_cell_glo,ind_cell_glo_glo)
546
         CALL bcast(ind_cell_glo_glo)
547
         ind_cell(:)=ind_cell_glo_glo(begin:end)
548
549
       ENDIF
550
       CALL Init_synchro_omp
551
552
!$OMP BARRIER
553
554
       IF (knon > 0) THEN
555
#ifdef CPP_VEGET
556
         CALL Init_intersurf(nbp_lon,nbp_lat,knon,ktindex,offset,orch_omp_size,orch_omp_rank,orch_comm,grid=grid_type)
557
#endif
558
       ENDIF
559
560
       CALL Synchro_omp
561
562
563
       IF (knon > 0) THEN
564
565
#ifdef CPP_VEGET
566
567
         CALL intersurf_initialize_gathered (itime+itau_phy-1, nbp_lon, nbp_lat, knon, ktindex, dtime, &
568
               lrestart_read, lrestart_write, lalo, contfrac, neighbours, resolution, date0, &
569
               zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, &
570
               cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
571
               precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
572
               evap, fluxsens, fluxlat, coastalflow, riverflow, &
573
               tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0m_new, &
574
               lon_scat, lat_scat, q2m(1:knon), t2m(1:knon), z0h_new(1:knon), nvm_orch, &
575
               grid=grid_type, bounds_latlon=bounds_lalo, cell_area=area, ind_cell_glo=ind_cell, &
576
               field_out_names=cfname_out, field_in_names=cfname_in(1:nbcf_in_orc), &
577
               coszang=yrmu0(1:knon))
578
#endif
579
       ENDIF
580
581
       CALL Synchro_omp
582
583
       albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
584
585
    ENDIF
586
587
!  swdown_vrai(1:knon) = swnet(1:knon)/(1. - albedo_keep(1:knon))
588
    swdown_vrai(1:knon) = swdown(1:knon)
589
!$OMP BARRIER
590
591
    IF (knon > 0) THEN
592
#ifdef CPP_VEGET
593
       IF (nvm_orch .NE. nvm_lmdz ) THEN
594
          abort_message='Pb de dimensiosn PFT: nvm_orch et nvm_lmdz differents.'
595
          CALL abort_physic(modname,abort_message,1)
596
       ENDIF
597
598
       CALL intersurf_main_gathered (itime+itau_phy, nbp_lon, nbp_lat, knon, ktindex, dtime,  &
599
            lrestart_read, lrestart_write, lalo, &
600
            contfrac, neighbours, resolution, date0, &
601
            zlev,  u1_lay(1:knon), v1_lay(1:knon), spechum(1:knon), temp_air(1:knon), epot_air(1:knon), ccanopy(1:knon), &
602
            cdrag(1:knon), petA_orc(1:knon), peqA_orc(1:knon), petB_orc(1:knon), peqB_orc(1:knon), &
603
            precip_rain(1:knon), precip_snow(1:knon), lwdown(1:knon), swnet(1:knon), swdown_vrai(1:knon), ps(1:knon), &
604
            evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
605
            tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0m_new(1:knon), &
606
            lon_scat, lat_scat, q2m(1:knon), t2m(1:knon), z0h_new(1:knon),&
607
            veget(1:knon,:),lai(1:knon,:),height(1:knon,:),&
608
            fields_out=yfields_out(1:knon,1:nbcf_out),  &
609
            fields_in=yfields_in(1:knon,1:nbcf_in_orc), &
610
            coszang=yrmu0(1:knon), run_off_lic=run_off_lic(1:knon), run_off_lic_frac=run_off_lic_frac(1:knon), blowingsnow_flux=blowingsnow_flux(1:knon))
611
#endif
612
    ENDIF
613
614
    CALL Synchro_omp
615
616
    albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
617
618
!* Send to coupler
619
!
620
    IF (type_ocean=='couple') THEN
621
       CALL cpl_send_land_fields(itime, knon, knindex, &
622
            riverflow, coastalflow)
623
       IF (landice_opt .GE. 2) THEN
624
          CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
625
       END IF
626
    ENDIF
627
628
    alb1_new(1:knon) = albedo_out(1:knon,1)
629
    alb2_new(1:knon) = albedo_out(1:knon,2)
630
631
! Convention orchidee: positif vers le haut
632
    fluxsens(1:knon) = -1. * fluxsens(1:knon)
633
    fluxlat(1:knon)  = -1. * fluxlat(1:knon)
634
635
!  evap     = -1. * evap
636
637
    IF (debut) lrestart_read = .FALSE.
638
639
    IF (debut) CALL Finalize_surf_para
640
641
! >> PC
642
! Decompressing variables into LMDz for the module carbon_cycle_mod
643
! nbcf_in can be zero, in which case the loop does not operate
644
! fields_in can then used elsewhere in the model
645
646
     fields_in(:,:)=0.0
647
648
     DO nb=1, nbcf_in_orc
649
       DO igrid = 1, knon
650
        ireal = knindex(igrid)
651
        fields_in(ireal,nb)=yfields_in(igrid,nb)
652
       ENDDO
653
       WRITE(*,*) 'surf_land_orchidee_mod --- yfields_in :',cfname_in(nb)
654
     ENDDO
655
! >> PC
656
657
  END SUBROUTINE surf_land_orchidee
658
!
659
!****************************************************************************************
660
!
661
  SUBROUTINE Init_orchidee_index(knon,knindex,offset,ktindex)
662
  USE mod_surf_para
663
  USE mod_grid_phy_lmdz
664
665
    INTEGER,INTENT(IN)    :: knon
666
    INTEGER,INTENT(IN)    :: knindex(klon)
667
    INTEGER,INTENT(OUT)   :: offset
668
    INTEGER,INTENT(OUT)   :: ktindex(klon)
669
670
    INTEGER               :: ktindex_glo(knon_glo)
671
    INTEGER               :: offset_para(0:omp_size*mpi_size-1)
672
    INTEGER               :: LastPoint
673
    INTEGER               :: task
674
675
    ktindex(1:knon)=knindex(1:knon)+(klon_mpi_begin-1)+(klon_omp_begin-1)+nbp_lon-1
676
677
    CALL gather_surf(ktindex(1:knon),ktindex_glo)
678
679
    IF (is_mpi_root .AND. is_omp_root) THEN
680
      LastPoint=0
681
      DO Task=0,mpi_size*omp_size-1
682
        IF (knon_glo_para(Task)>0) THEN
683
           offset_para(task)= LastPoint-MOD(LastPoint,nbp_lon)
684
           LastPoint=ktindex_glo(knon_glo_end_para(task))
685
        ENDIF
686
      ENDDO
687
    ENDIF
688
689
    CALL bcast(offset_para)
690
691
    offset=offset_para(omp_size*mpi_rank+omp_rank)
692
693
    ktindex(1:knon)=ktindex(1:knon)-offset
694
695
  END SUBROUTINE Init_orchidee_index
696
697
!
698
!************************* ***************************************************************
699
!
700
701
  SUBROUTINE Get_orchidee_communicator(orch_comm, orch_mpi_size, orch_mpi_rank, orch_omp_size,orch_omp_rank)
702
  USE  mod_surf_para
703
704
#ifdef CPP_MPI
705
    INCLUDE 'mpif.h'
706
#endif
707
708
    INTEGER,INTENT(OUT) :: orch_comm
709
    INTEGER,INTENT(OUT) :: orch_mpi_size
710
    INTEGER,INTENT(OUT) :: orch_mpi_rank
711
    INTEGER,INTENT(OUT) :: orch_omp_size
712
    INTEGER,INTENT(OUT) :: orch_omp_rank
713
    INTEGER             :: color
714
    INTEGER             :: i,ierr
715
!
716
! End definition
717
!****************************************************************************************
718
719
    IF (is_omp_root) THEN
720
721
      IF (knon_mpi==0) THEN
722
         color = 0
723
      ELSE
724
         color = 1
725
      ENDIF
726
727
#ifdef CPP_MPI
728
      CALL MPI_COMM_SPLIT(COMM_LMDZ_PHY,color,mpi_rank,orch_comm,ierr)
729
      CALL MPI_COMM_SIZE(orch_comm,orch_mpi_size,ierr)
730
      CALL MPI_COMM_RANK(orch_comm,orch_mpi_rank,ierr)
731
#endif
732
733
    ENDIF
734
    CALL bcast_omp(orch_comm)
735
736
    IF (knon_mpi /= 0) THEN
737
      orch_omp_size=0
738
      DO i=0,omp_size-1
739
        IF (knon_omp_para(i) /=0) THEN
740
          orch_omp_size=orch_omp_size+1
741
          IF (i==omp_rank) orch_omp_rank=orch_omp_size-1
742
        ENDIF
743
      ENDDO
744
    ENDIF
745
746
  END SUBROUTINE Get_orchidee_communicator
747
!
748
!****************************************************************************************
749
!
750
751
  SUBROUTINE Init_neighbours(knon,neighbours,knindex,pctsrf)
752
    USE mod_grid_phy_lmdz
753
    USE mod_surf_para
754
    USE indice_sol_mod
755
756
#ifdef CPP_MPI
757
    INCLUDE 'mpif.h'
758
#endif
759
760
! Input arguments
761
!****************************************************************************************
762
    INTEGER, INTENT(IN)                     :: knon
763
    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
764
    REAL, DIMENSION(klon), INTENT(IN)       :: pctsrf
765
766
! Output arguments
767
!****************************************************************************************
768
    INTEGER, DIMENSION(knon,8), INTENT(OUT) :: neighbours
769
770
! Local variables
771
!****************************************************************************************
772
    INTEGER                              :: i, igrid, jj, ij, iglob
773
    INTEGER                              :: ierr, ireal, index
774
    INTEGER, DIMENSION(8,3)              :: off_ini
775
    INTEGER, DIMENSION(8)                :: offset
776
    INTEGER, DIMENSION(nbp_lon,nbp_lat)  :: correspond
777
    INTEGER, DIMENSION(knon_glo)         :: ktindex_glo
778
    INTEGER, DIMENSION(knon_glo,8)       :: neighbours_glo
779
    REAL, DIMENSION(klon_glo)            :: pctsrf_glo
780
    INTEGER                              :: ktindex(klon)
781
!
782
! End definition
783
!****************************************************************************************
784
785
    ktindex(1:knon)=knindex(1:knon)+(klon_mpi_begin-1)+(klon_omp_begin-1)+nbp_lon-1
786
787
    CALL gather_surf(ktindex(1:knon),ktindex_glo)
788
    CALL gather(pctsrf,pctsrf_glo)
789
790
    IF (is_mpi_root .AND. is_omp_root) THEN
791
      neighbours_glo(:,:)=-1
792
!  Initialisation des offset
793
!
794
! offset bord ouest
795
       off_ini(1,1) = - nbp_lon   ; off_ini(2,1) = - nbp_lon + 1     ; off_ini(3,1) = 1
796
       off_ini(4,1) = nbp_lon + 1 ; off_ini(5,1) = nbp_lon           ; off_ini(6,1) = 2 * nbp_lon - 1
797
       off_ini(7,1) = nbp_lon -1  ; off_ini(8,1) = - 1
798
! offset point normal
799
       off_ini(1,2) = - nbp_lon   ; off_ini(2,2) = - nbp_lon + 1     ; off_ini(3,2) = 1
800
       off_ini(4,2) = nbp_lon + 1 ; off_ini(5,2) = nbp_lon           ; off_ini(6,2) = nbp_lon - 1
801
       off_ini(7,2) = -1          ; off_ini(8,2) = - nbp_lon - 1
802
! offset bord   est
803
       off_ini(1,3) = - nbp_lon   ; off_ini(2,3) = - 2 * nbp_lon + 1 ; off_ini(3,3) = - nbp_lon + 1
804
       off_ini(4,3) =  1          ; off_ini(5,3) = nbp_lon           ; off_ini(6,3) = nbp_lon - 1
805
       off_ini(7,3) = -1          ; off_ini(8,3) = - nbp_lon - 1
806
!
807
! Attention aux poles
808
!
809
       DO igrid = 1, knon_glo
810
          index = ktindex_glo(igrid)
811
          jj = INT((index - 1)/nbp_lon) + 1
812
          ij = index - (jj - 1) * nbp_lon
813
          correspond(ij,jj) = igrid
814
       ENDDO
815
!sonia : Les mailles des voisines doivent etre toutes egales (pour couplage orchidee)
816
       IF (knon_glo == 1) THEN
817
         igrid = 1
818
         DO i = 1,8
819
           neighbours_glo(igrid, i) = igrid
820
         ENDDO
821
       ELSE
822
823
       DO igrid = 1, knon_glo
824
          iglob = ktindex_glo(igrid)
825
826
          IF (MOD(iglob, nbp_lon) == 1) THEN
827
             offset = off_ini(:,1)
828
          ELSE IF(MOD(iglob, nbp_lon) == 0) THEN
829
             offset = off_ini(:,3)
830
          ELSE
831
             offset = off_ini(:,2)
832
          ENDIF
833
834
          DO i = 1, 8
835
             index = iglob + offset(i)
836
             ireal = (MIN(MAX(1, index - nbp_lon + 1), klon_glo))
837
             IF (pctsrf_glo(ireal) > EPSFRA) THEN
838
                jj = INT((index - 1)/nbp_lon) + 1
839
                ij = index - (jj - 1) * nbp_lon
840
                neighbours_glo(igrid, i) = correspond(ij, jj)
841
             ENDIF
842
          ENDDO
843
       ENDDO
844
       ENDIF !fin knon_glo == 1
845
846
    ENDIF
847
848
    DO i = 1, 8
849
      CALL scatter_surf(neighbours_glo(:,i),neighbours(1:knon,i))
850
    ENDDO
851
  END SUBROUTINE Init_neighbours
852
853
!
854
!****************************************************************************************
855
!
856
#endif
857
#endif
858
#endif
859
#endif
860
#endif
861
END MODULE surf_land_orchidee_mod