source: LMDZ5/trunk/libf/phylmd/surf_land_orchidee_mod.F90 @ 2351

Last change on this file since 2351 was 2351, checked in by Ehouarn Millour, 9 years ago

More on physics/dynamics separation and cleanup:

  • Set things up so that all physics-related initializations are done via iniphysiq.
  • Created a "geometry_mod.F90" module in phy_common to store information on the loacl grid (i.e. replaces comgeomphy) and moreover give these variables more obvious names (e.g.: rlond => longitude, rlatd => latitude, airephy => cell_area).
  • removed obsolete comgeomphy.h and comgeomphy.F90

EM

  • 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: 22.2 KB
Line 
1!
2MODULE surf_land_orchidee_mod
3#ifndef ORCHIDEE_NOOPENMP
4!
5! This module controles the interface towards the model ORCHIDEE
6!
7! Subroutines in this module : surf_land_orchidee
8!                              Init_orchidee_index
9!                              Get_orchidee_communicator
10!                              Init_neighbours
11
12  USE dimphy
13#ifdef CPP_VEGET
14  USE intersurf     ! module d'ORCHIDEE
15#endif
16  USE cpl_mod,      ONLY : cpl_send_land_fields
17  USE surface_data, ONLY : type_ocean
18  USE geometry_mod, ONLY : dx, dy
19  USE mod_grid_phy_lmdz
20  USE mod_phys_lmdz_para, mpi_root_rank=>mpi_master
21
22  IMPLICIT NONE
23
24  PRIVATE
25  PUBLIC  :: surf_land_orchidee
26
27  LOGICAL, ALLOCATABLE, SAVE :: flag_omp(:)
28CONTAINS
29!
30!****************************************************************************************
31
32  SUBROUTINE surf_land_orchidee(itime, dtime, date0, knon, &
33       knindex, rlon, rlat, pctsrf, &
34       debut, lafin, &
35       plev,  u1_lay, v1_lay, gustiness, temp_air, spechum, epot_air, ccanopy, &
36       tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
37       precip_rain, precip_snow, lwdown, swnet, swdown, &
38       ps, q2m, t2m, &
39       evap, fluxsens, fluxlat, &             
40       tsol_rad, tsurf_new, alb1_new, alb2_new, &
41       emis_new, z0_new, qsurf)
42
43    USE mod_surf_para
44    USE mod_synchro_omp
45    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl
46    USE indice_sol_mod
47    USE print_control_mod, ONLY: lunout
48    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
49#ifdef CPP_VEGET
50    USE time_phylmdz_mod, ONLY: itau_phy
51#endif
52!   
53! Cette routine sert d'interface entre le modele atmospherique et le
54! modele de sol continental. Appel a sechiba
55!
56! L. Fairhead 02/2000
57!
58! input:
59!   itime        numero du pas de temps
60!   dtime        pas de temps de la physique (en s)
61!   nisurf       index de la surface a traiter (1 = sol continental)
62!   knon         nombre de points de la surface a traiter
63!   knindex      index des points de la surface a traiter
64!   rlon         longitudes de la grille entiere
65!   rlat         latitudes de la grille entiere
66!   pctsrf       tableau des fractions de surface de chaque maille
67!   debut        logical: 1er appel a la physique (lire les restart)
68!   lafin        logical: dernier appel a la physique (ecrire les restart)
69!                     (si false calcul simplifie des fluxs sur les continents)
70!   plev         hauteur de la premiere couche (Pa)     
71!   u1_lay       vitesse u 1ere couche
72!   v1_lay       vitesse v 1ere couche
73!   temp_air     temperature de l'air 1ere couche
74!   spechum      humidite specifique 1ere couche
75!   epot_air     temp pot de l'air
76!   ccanopy      concentration CO2 canopee, correspond au co2_send de
77!                carbon_cycle_mod ou valeur constant co2_ppm
78!   tq_cdrag     cdrag
79!   petAcoef     coeff. A de la resolution de la CL pour t
80!   peqAcoef     coeff. A de la resolution de la CL pour q
81!   petBcoef     coeff. B de la resolution de la CL pour t
82!   peqBcoef     coeff. B de la resolution de la CL pour q
83!   precip_rain  precipitation liquide
84!   precip_snow  precipitation solide
85!   lwdown       flux IR descendant a la surface
86!   swnet        flux solaire net
87!   swdown       flux solaire entrant a la surface
88!   ps           pression au sol
89!   radsol       rayonnement net aus sol (LW + SW)
90!   
91!
92! output:
93!   evap         evaporation totale
94!   fluxsens     flux de chaleur sensible
95!   fluxlat      flux de chaleur latente
96!   tsol_rad     
97!   tsurf_new    temperature au sol
98!   alb1_new     albedo in visible SW interval
99!   alb2_new     albedo in near IR interval
100!   emis_new     emissivite
101!   z0_new       surface roughness
102!   qsurf        air moisture at surface
103!
104    INCLUDE "YOMCST.h"
105 
106!
107! Parametres d'entree
108!****************************************************************************************
109    INTEGER, INTENT(IN)                       :: itime
110    REAL, INTENT(IN)                          :: dtime
111    REAL, INTENT(IN)                          :: date0
112    INTEGER, INTENT(IN)                       :: knon
113    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex
114    LOGICAL, INTENT(IN)                       :: debut, lafin
115    REAL, DIMENSION(klon,nbsrf), INTENT(IN)   :: pctsrf
116    REAL, DIMENSION(klon), INTENT(IN)         :: rlon, rlat
117    REAL, DIMENSION(klon), INTENT(IN)         :: plev
118    REAL, DIMENSION(klon), INTENT(IN)         :: u1_lay, v1_lay, gustiness
119    REAL, DIMENSION(klon), INTENT(IN)         :: temp_air, spechum
120    REAL, DIMENSION(klon), INTENT(IN)         :: epot_air, ccanopy
121    REAL, DIMENSION(klon), INTENT(IN)         :: tq_cdrag
122    REAL, DIMENSION(klon), INTENT(IN)         :: petAcoef, peqAcoef
123    REAL, DIMENSION(klon), INTENT(IN)         :: petBcoef, peqBcoef
124    REAL, DIMENSION(klon), INTENT(IN)         :: precip_rain, precip_snow
125    REAL, DIMENSION(klon), INTENT(IN)         :: lwdown, swnet, swdown, ps
126    REAL, DIMENSION(klon), INTENT(IN)         :: q2m, t2m
127
128! Parametres de sortie
129!****************************************************************************************
130    REAL, DIMENSION(klon), INTENT(OUT)        :: evap, fluxsens, fluxlat, qsurf
131    REAL, DIMENSION(klon), INTENT(OUT)        :: tsol_rad, tsurf_new
132    REAL, DIMENSION(klon), INTENT(OUT)        :: alb1_new, alb2_new
133    REAL, DIMENSION(klon), INTENT(OUT)        :: emis_new, z0_new
134
135! Local
136!****************************************************************************************
137    INTEGER                                   :: ij, jj, igrid, ireal, index
138    INTEGER                                   :: error
139    REAL, DIMENSION(klon)                     :: swdown_vrai
140    CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
141    CHARACTER (len = 80)                      :: abort_message
142    LOGICAL,SAVE                              :: check = .FALSE.
143    !$OMP THREADPRIVATE(check)
144
145! type de couplage dans sechiba
146!  character (len=10)   :: coupling = 'implicit'
147! drapeaux controlant les appels dans SECHIBA
148!  type(control_type), save   :: control_in
149! Preserved albedo
150    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: albedo_keep, zlev
151    !$OMP THREADPRIVATE(albedo_keep,zlev)
152! coordonnees geographiques
153    REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: lalo
154    !$OMP THREADPRIVATE(lalo)
155! pts voisins
156    INTEGER,ALLOCATABLE, DIMENSION(:,:), SAVE :: neighbours
157    !$OMP THREADPRIVATE(neighbours)
158! fractions continents
159    REAL,ALLOCATABLE, DIMENSION(:), SAVE      :: contfrac
160    !$OMP THREADPRIVATE(contfrac)
161! resolution de la grille
162    REAL, ALLOCATABLE, DIMENSION (:,:), SAVE  :: resolution
163    !$OMP THREADPRIVATE(resolution)
164
165    REAL, ALLOCATABLE, DIMENSION (:,:), SAVE  :: lon_scat, lat_scat 
166    !$OMP THREADPRIVATE(lon_scat,lat_scat)
167
168    LOGICAL, SAVE                             :: lrestart_read = .TRUE.
169    !$OMP THREADPRIVATE(lrestart_read)
170    LOGICAL, SAVE                             :: lrestart_write = .FALSE.
171    !$OMP THREADPRIVATE(lrestart_write)
172
173    REAL, DIMENSION(knon,2)                   :: albedo_out
174
175! Pb de nomenclature
176    REAL, DIMENSION(klon)                     :: petA_orc, peqA_orc
177    REAL, DIMENSION(klon)                     :: petB_orc, peqB_orc
178! Pb de correspondances de grilles
179    INTEGER, DIMENSION(:), SAVE, ALLOCATABLE  :: ig, jg
180    !$OMP THREADPRIVATE(ig,jg)
181    INTEGER :: indi, indj
182    INTEGER, SAVE, ALLOCATABLE,DIMENSION(:)   :: ktindex
183    !$OMP THREADPRIVATE(ktindex)
184
185! Essai cdrag
186    REAL, DIMENSION(klon)                     :: cdrag
187    INTEGER,SAVE                              :: offset
188    !$OMP THREADPRIVATE(offset)
189
190    REAL, DIMENSION(klon_glo)                 :: rlon_g,rlat_g
191    INTEGER, SAVE                             :: orch_comm
192    !$OMP THREADPRIVATE(orch_comm)
193
194    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: coastalflow
195    !$OMP THREADPRIVATE(coastalflow)
196    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: riverflow
197    !$OMP THREADPRIVATE(riverflow)
198   
199    INTEGER :: orch_omp_rank
200    INTEGER :: orch_omp_size
201!
202! Fin definition
203!****************************************************************************************
204
205    IF (check) WRITE(lunout,*)'Entree ', modname
206 
207! Initialisation
208 
209    IF (debut) THEN
210! Test of coherence between variable ok_veget and cpp key CPP_VEGET
211#ifndef CPP_VEGET
212       abort_message='Pb de coherence: ok_veget = .true. mais CPP_VEGET = .false.'
213       CALL abort_physic(modname,abort_message,1)
214#endif
215
216       CALL Init_surf_para(knon)
217       ALLOCATE(ktindex(knon))
218       IF ( .NOT. ALLOCATED(albedo_keep)) THEN
219!ym          ALLOCATE(albedo_keep(klon))
220!ym bizarre que non alloué en knon precedement
221          ALLOCATE(albedo_keep(knon))
222          ALLOCATE(zlev(knon))
223       ENDIF
224! Pb de correspondances de grilles
225       ALLOCATE(ig(klon))
226       ALLOCATE(jg(klon))
227       ig(1) = 1
228       jg(1) = 1
229       indi = 0
230       indj = 2
231       DO igrid = 2, klon - 1
232          indi = indi + 1
233          IF ( indi > nbp_lon) THEN
234             indi = 1
235             indj = indj + 1
236          ENDIF
237          ig(igrid) = indi
238          jg(igrid) = indj
239       ENDDO
240       ig(klon) = 1
241       jg(klon) = nbp_lat
242
243       IF ((.NOT. ALLOCATED(lalo))) THEN
244          ALLOCATE(lalo(knon,2), stat = error)
245          IF (error /= 0) THEN
246             abort_message='Pb allocation lalo'
247             CALL abort_physic(modname,abort_message,1)
248          ENDIF
249       ENDIF
250       IF ((.NOT. ALLOCATED(lon_scat))) THEN
251          ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
252          IF (error /= 0) THEN
253             abort_message='Pb allocation lon_scat'
254             CALL abort_physic(modname,abort_message,1)
255          ENDIF
256       ENDIF
257       IF ((.NOT. ALLOCATED(lat_scat))) THEN
258          ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
259          IF (error /= 0) THEN
260             abort_message='Pb allocation lat_scat'
261             CALL abort_physic(modname,abort_message,1)
262          ENDIF
263       ENDIF
264       lon_scat = 0.
265       lat_scat = 0.
266       DO igrid = 1, knon
267          index = knindex(igrid)
268          lalo(igrid,2) = rlon(index)
269          lalo(igrid,1) = rlat(index)
270       ENDDO
271
272       
273       
274       CALL Gather(rlon,rlon_g)
275       CALL Gather(rlat,rlat_g)
276
277       IF (is_mpi_root) THEN
278          index = 1
279          DO jj = 2, nbp_lat-1
280             DO ij = 1, nbp_lon
281                index = index + 1
282                lon_scat(ij,jj) = rlon_g(index)
283                lat_scat(ij,jj) = rlat_g(index)
284             ENDDO
285          ENDDO
286          lon_scat(:,1) = lon_scat(:,2)
287          lat_scat(:,1) = rlat_g(1)
288          lon_scat(:,nbp_lat) = lon_scat(:,2)
289          lat_scat(:,nbp_lat) = rlat_g(klon_glo)
290       ENDIF
291   
292       CALL bcast(lon_scat)
293       CALL bcast(lat_scat)
294!
295! Allouer et initialiser le tableau des voisins et des fraction de continents
296!
297       IF ( (.NOT.ALLOCATED(neighbours))) THEN
298          ALLOCATE(neighbours(knon,8), stat = error)
299          IF (error /= 0) THEN
300             abort_message='Pb allocation neighbours'
301             CALL abort_physic(modname,abort_message,1)
302          ENDIF
303       ENDIF
304       neighbours = -1.
305       IF (( .NOT. ALLOCATED(contfrac))) THEN
306          ALLOCATE(contfrac(knon), stat = error)
307          IF (error /= 0) THEN
308             abort_message='Pb allocation contfrac'
309             CALL abort_physic(modname,abort_message,1)
310          ENDIF
311       ENDIF
312
313       DO igrid = 1, knon
314          ireal = knindex(igrid)
315          contfrac(igrid) = pctsrf(ireal,is_ter)
316       ENDDO
317
318
319       CALL Init_neighbours(knon,neighbours,knindex,pctsrf(:,is_ter))
320
321!
322!  Allocation et calcul resolutions
323       IF ( (.NOT.ALLOCATED(resolution))) THEN
324          ALLOCATE(resolution(knon,2), stat = error)
325          IF (error /= 0) THEN
326             abort_message='Pb allocation resolution'
327             CALL abort_physic(modname,abort_message,1)
328          ENDIF
329       ENDIF
330       DO igrid = 1, knon
331          ij = knindex(igrid)
332          resolution(igrid,1) = dx(ij)
333          resolution(igrid,2) = dy(ij)
334       ENDDO
335     
336       ALLOCATE(coastalflow(klon), stat = error)
337       IF (error /= 0) THEN
338          abort_message='Pb allocation coastalflow'
339          CALL abort_physic(modname,abort_message,1)
340       ENDIF
341       
342       ALLOCATE(riverflow(klon), stat = error)
343       IF (error /= 0) THEN
344          abort_message='Pb allocation riverflow'
345          CALL abort_physic(modname,abort_message,1)
346       ENDIF
347!
348! carbon_cycle_cpl not possible with this interface and version of ORHCHIDEE
349!
350       IF (carbon_cycle_cpl) THEN
351          abort_message='carbon_cycle_cpl not yet possible with this interface of ORCHIDEE'
352          CALL abort_physic(modname,abort_message,1)
353       END IF
354       
355    ENDIF                          ! (fin debut)
356 
357
358!
359! Appel a la routine sols continentaux
360!
361    IF (lafin) lrestart_write = .TRUE.
362    IF (check) WRITE(lunout,*)'lafin ',lafin,lrestart_write
363     
364    petA_orc(1:knon) = petBcoef(1:knon) * dtime
365    petB_orc(1:knon) = petAcoef(1:knon)
366    peqA_orc(1:knon) = peqBcoef(1:knon) * dtime
367    peqB_orc(1:knon) = peqAcoef(1:knon)
368
369    cdrag = 0.
370    cdrag(1:knon) = tq_cdrag(1:knon)
371
372! zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/287.05*temp_air(1:knon))*9.80665)
373!    zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/RD*temp_air(1:knon))*RG)
374     zlev(1:knon) = plev(1:knon)*RD*temp_air(1:knon)/((ps(1:knon)*100.0)*RG)
375
376
377! PF et PASB
378!   where(cdrag > 0.01)
379!     cdrag = 0.01
380!   endwhere
381!  write(*,*)'Cdrag = ',minval(cdrag),maxval(cdrag)
382
383 
384    IF (debut) THEN
385       CALL Init_orchidee_index(knon,knindex,offset,ktindex)
386       CALL Get_orchidee_communicator(orch_comm,orch_omp_size,orch_omp_rank)
387       CALL Init_synchro_omp
388       
389       IF (knon > 0) THEN
390#ifdef CPP_VEGET
391         CALL Init_intersurf(nbp_lon,nbp_lat,knon,ktindex,offset,orch_omp_size,orch_omp_rank,orch_comm)
392#endif
393       ENDIF
394
395       
396       IF (knon > 0) THEN
397
398#ifdef CPP_VEGET
399          CALL intersurf_main (itime+itau_phy-1, nbp_lon, nbp_lat, knon, ktindex, dtime, &
400               lrestart_read, lrestart_write, lalo, &
401               contfrac, neighbours, resolution, date0, &
402               zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
403               cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
404               precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
405               evap, fluxsens, fluxlat, coastalflow, riverflow, &
406               tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
407               lon_scat, lat_scat, q2m, t2m)
408#endif         
409       ENDIF
410
411       CALL Synchro_omp
412
413       albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
414
415    ENDIF
416
417   
418!  swdown_vrai(1:knon) = swnet(1:knon)/(1. - albedo_keep(1:knon))
419    swdown_vrai(1:knon) = swdown(1:knon)
420
421    IF (knon > 0) THEN
422#ifdef CPP_VEGET   
423       CALL intersurf_main (itime+itau_phy, nbp_lon, nbp_lat, knon, ktindex, dtime,  &
424            lrestart_read, lrestart_write, lalo, &
425            contfrac, neighbours, resolution, date0, &
426            zlev,  u1_lay(1:knon), v1_lay(1:knon), spechum(1:knon), temp_air(1:knon), epot_air(1:knon), ccanopy(1:knon), &
427            cdrag(1:knon), petA_orc(1:knon), peqA_orc(1:knon), petB_orc(1:knon), peqB_orc(1:knon), &
428            precip_rain(1:knon), precip_snow(1:knon), lwdown(1:knon), swnet(1:knon), swdown_vrai(1:knon), ps(1:knon), &
429            evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
430            tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
431            lon_scat, lat_scat, q2m, t2m)
432#endif       
433    ENDIF
434
435    CALL Synchro_omp
436   
437    albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
438
439!* Send to coupler
440!
441    IF (type_ocean=='couple') THEN
442       CALL cpl_send_land_fields(itime, knon, knindex, &
443            riverflow, coastalflow)
444    ENDIF
445
446    alb1_new(1:knon) = albedo_out(1:knon,1)
447    alb2_new(1:knon) = albedo_out(1:knon,2)
448
449! Convention orchidee: positif vers le haut
450    fluxsens(1:knon) = -1. * fluxsens(1:knon)
451    fluxlat(1:knon)  = -1. * fluxlat(1:knon)
452   
453!  evap     = -1. * evap
454
455    IF (debut) lrestart_read = .FALSE.
456   
457    IF (debut) CALL Finalize_surf_para
458
459   
460  END SUBROUTINE surf_land_orchidee
461!
462!****************************************************************************************
463!
464  SUBROUTINE Init_orchidee_index(knon,knindex,offset,ktindex)
465  USE mod_surf_para
466  USE mod_grid_phy_lmdz
467 
468    INTEGER,INTENT(IN)    :: knon
469    INTEGER,INTENT(IN)    :: knindex(klon)   
470    INTEGER,INTENT(OUT)   :: offset
471    INTEGER,INTENT(OUT)   :: ktindex(klon)
472   
473    INTEGER               :: ktindex_glo(knon_glo)
474    INTEGER               :: offset_para(0:omp_size*mpi_size-1)
475    INTEGER               :: LastPoint
476    INTEGER               :: task
477   
478    ktindex(1:knon)=knindex(1:knon)+(klon_mpi_begin-1)+(klon_omp_begin-1)+nbp_lon-1
479   
480    CALL gather_surf(ktindex(1:knon),ktindex_glo)
481   
482    IF (is_mpi_root .AND. is_omp_root) THEN
483      LastPoint=0
484      DO Task=0,mpi_size*omp_size-1
485        IF (knon_glo_para(Task)>0) THEN
486           offset_para(task)= LastPoint-MOD(LastPoint,nbp_lon)
487           LastPoint=ktindex_glo(knon_glo_end_para(task))
488        ENDIF
489      ENDDO
490    ENDIF
491   
492    CALL bcast(offset_para)
493   
494    offset=offset_para(omp_size*mpi_rank+omp_rank)
495   
496    ktindex(1:knon)=ktindex(1:knon)-offset
497
498  END SUBROUTINE Init_orchidee_index
499
500!
501!************************* ***************************************************************
502!
503
504  SUBROUTINE Get_orchidee_communicator(orch_comm,orch_omp_size,orch_omp_rank)
505  USE  mod_surf_para
506     
507#ifdef CPP_MPI
508    INCLUDE 'mpif.h'
509#endif   
510
511    INTEGER,INTENT(OUT) :: orch_comm
512    INTEGER,INTENT(OUT) :: orch_omp_size
513    INTEGER,INTENT(OUT) :: orch_omp_rank
514    INTEGER             :: color
515    INTEGER             :: i,ierr
516!
517! End definition
518!****************************************************************************************
519   
520   
521    IF (is_omp_root) THEN         
522     
523      IF (knon_mpi==0) THEN
524         color = 0
525      ELSE
526         color = 1
527      ENDIF
528   
529#ifdef CPP_MPI   
530      CALL MPI_COMM_SPLIT(COMM_LMDZ_PHY,color,mpi_rank,orch_comm,ierr)
531#endif
532   
533    ENDIF
534    CALL bcast_omp(orch_comm)
535   
536    IF (knon_mpi /= 0) THEN
537      orch_omp_size=0
538      DO i=0,omp_size-1
539        IF (knon_omp_para(i) /=0) THEN
540          orch_omp_size=orch_omp_size+1
541          IF (i==omp_rank) orch_omp_rank=orch_omp_size-1
542        ENDIF
543      ENDDO
544    ENDIF
545   
546   
547  END SUBROUTINE Get_orchidee_communicator
548!
549!****************************************************************************************
550
551
552  SUBROUTINE Init_neighbours(knon,neighbours,knindex,pctsrf)
553    USE mod_grid_phy_lmdz
554    USE mod_surf_para   
555    USE indice_sol_mod
556
557#ifdef CPP_MPI
558    INCLUDE 'mpif.h'
559#endif   
560
561! Input arguments
562!****************************************************************************************
563    INTEGER, INTENT(IN)                     :: knon
564    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
565    REAL, DIMENSION(klon), INTENT(IN)       :: pctsrf
566   
567! Output arguments
568!****************************************************************************************
569    INTEGER, DIMENSION(knon,8), INTENT(OUT) :: neighbours
570
571! Local variables
572!****************************************************************************************
573    INTEGER                              :: i, igrid, jj, ij, iglob
574    INTEGER                              :: ierr, ireal, index
575    INTEGER, DIMENSION(8,3)              :: off_ini
576    INTEGER, DIMENSION(8)                :: offset 
577    INTEGER, DIMENSION(nbp_lon,nbp_lat)  :: correspond
578    INTEGER, DIMENSION(knon_glo)         :: ktindex_glo
579    INTEGER, DIMENSION(knon_glo,8)       :: neighbours_glo
580    REAL, DIMENSION(klon_glo)            :: pctsrf_glo
581    INTEGER                              :: ktindex(klon)
582!
583! End definition
584!****************************************************************************************
585
586    ktindex(1:knon)=knindex(1:knon)+(klon_mpi_begin-1)+(klon_omp_begin-1)+nbp_lon-1
587   
588    CALL gather_surf(ktindex(1:knon),ktindex_glo)
589    CALL gather(pctsrf,pctsrf_glo)
590   
591    IF (is_mpi_root .AND. is_omp_root) THEN
592      neighbours_glo(:,:)=-1
593!  Initialisation des offset   
594!
595! offset bord ouest
596       off_ini(1,1) = - nbp_lon   ; off_ini(2,1) = - nbp_lon + 1     ; off_ini(3,1) = 1
597       off_ini(4,1) = nbp_lon + 1 ; off_ini(5,1) = nbp_lon           ; off_ini(6,1) = 2 * nbp_lon - 1
598       off_ini(7,1) = nbp_lon -1  ; off_ini(8,1) = - 1
599! offset point normal
600       off_ini(1,2) = - nbp_lon   ; off_ini(2,2) = - nbp_lon + 1     ; off_ini(3,2) = 1
601       off_ini(4,2) = nbp_lon + 1 ; off_ini(5,2) = nbp_lon           ; off_ini(6,2) = nbp_lon - 1
602       off_ini(7,2) = -1          ; off_ini(8,2) = - nbp_lon - 1
603! offset bord   est
604       off_ini(1,3) = - nbp_lon   ; off_ini(2,3) = - 2 * nbp_lon + 1 ; off_ini(3,3) = - nbp_lon + 1
605       off_ini(4,3) =  1          ; off_ini(5,3) = nbp_lon           ; off_ini(6,3) = nbp_lon - 1
606       off_ini(7,3) = -1          ; off_ini(8,3) = - nbp_lon - 1
607!
608!
609! Attention aux poles
610!
611       DO igrid = 1, knon_glo
612          index = ktindex_glo(igrid)
613          jj = INT((index - 1)/nbp_lon) + 1
614          ij = index - (jj - 1) * nbp_lon
615          correspond(ij,jj) = igrid
616       ENDDO
617!sonia : Les mailles des voisines doivent etre toutes egales (pour couplage orchidee)
618       IF (knon_glo == 1) THEN
619         igrid = 1
620         DO i = 1,8
621           neighbours_glo(igrid, i) = igrid
622         ENDDO
623       ELSE
624       print*,'sonia : knon_glo,ij,jj', knon_glo, ij,jj
625       
626       DO igrid = 1, knon_glo
627          iglob = ktindex_glo(igrid)
628         
629          IF (MOD(iglob, nbp_lon) == 1) THEN
630             offset = off_ini(:,1)
631          ELSE IF(MOD(iglob, nbp_lon) == 0) THEN
632             offset = off_ini(:,3)
633          ELSE
634             offset = off_ini(:,2)
635          ENDIF
636         
637          DO i = 1, 8
638             index = iglob + offset(i)
639             ireal = (MIN(MAX(1, index - nbp_lon + 1), klon_glo))
640             IF (pctsrf_glo(ireal) > EPSFRA) THEN
641                jj = INT((index - 1)/nbp_lon) + 1
642                ij = index - (jj - 1) * nbp_lon
643                neighbours_glo(igrid, i) = correspond(ij, jj)
644             ENDIF
645          ENDDO
646       ENDDO
647       ENDIF !fin knon_glo == 1
648
649    ENDIF
650   
651    DO i = 1, 8
652      CALL scatter_surf(neighbours_glo(:,i),neighbours(1:knon,i))
653    ENDDO
654  END SUBROUTINE Init_neighbours
655
656!
657!****************************************************************************************
658!
659#endif
660END MODULE surf_land_orchidee_mod
Note: See TracBrowser for help on using the repository browser.