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

Last change on this file since 2271 was 2240, checked in by fhourdin, 10 years ago

Revisite de la formule des flux de surface
(en priorité sur l'océan) en tenant compte des bourrasques de
vent et de la différence entre les hauteurs de rugosités pour
la quantité de mouvement, l'enthalpie et éventuellement l'humidité.

Etape 1 :
Introduction d'un calcul de gustiness dans la physique
gustiness(:)=f_gust_bl * ale_bl + f_gust_wk * ame_wk
Cette variable est passée ensuite jusqu'au fin fond de la couche limite.
L'étape 1 est prête à commettre, ne nécessite pas de nouvelles
variables dans les startphy et assure la convergence numérique.

Introduction of gustiness in the surface flux computation.
Gustiness is computed from as
gustiness(:)=f_gust_bl * ale_bl + f_gust_wk * ame_wk
and pass through pbl_surface down to the routines that compute
surface fluxes.

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