source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/surf_land_orchidee_mod.F90

Last change on this file was 3331, checked in by acozic, 6 years ago

Add modification for isotopes

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