source: LMDZ6/trunk/libf/phylmd/surf_land_orchidee_mod.F90 @ 5433

Last change on this file since 5433 was 5305, checked in by abarral, 2 months ago

Turn YOMCST2.h.h into module

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