Changeset 3413


Ignore:
Timestamp:
Nov 12, 2018, 1:52:29 PM (6 years ago)
Author:
Laurent Fairhead
Message:

Inclusion of Yann's latest (summer/fall 2018) modifications for
convergence of DYNAMICO/LMDZ physics
YM/LF

Location:
LMDZ6/branches/DYNAMICO-conv
Files:
23 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/DYNAMICO-conv/DefLists/context_lmdz.xml

    r3411 r3413  
    1212  <file_definition src="./file_def_histhf_lmdz.xml"/>
    1313  <file_definition src="./file_def_histins_lmdz.xml"/>
    14 <!--  <file_definition src="./file_def_histLES_lmdz.xml"/> -->
     14  <file_definition src="./file_def_histLES_lmdz.xml"/>
    1515  <file_definition src="./file_def_histmth_lmdz.xml"/>
    16 <!--  <file_definition src="./file_def_histstn_lmdz.xml"/> -->
     16  <file_definition src="./file_def_histstn_lmdz.xml"/>
    1717  <file_definition src="./file_def_histmthNMC_lmdz.xml"/>
    1818  <file_definition src="./file_def_histdayNMC_lmdz.xml"/>
    1919  <file_definition src="./file_def_histhfNMC_lmdz.xml"/>
    20 <!--  <file_definition src="./file_def_histmthCOSP_lmdz.xml"/> -->
    21 <!--  <file_definition src="./file_def_histdayCOSP_lmdz.xml"/> -->
    22 <!--  <file_definition src="./file_def_histhfCOSP_lmdz.xml"/> -->
     20  <file_definition src="./file_def_histmthCOSP_lmdz.xml"/> -->
     21  <file_definition src="./file_def_histdayCOSP_lmdz.xml"/> -->
     22  <file_definition src="./file_def_histhfCOSP_lmdz.xml"/> -->
    2323  <file_definition src="./file_def_histstrataer_lmdz.xml"/>
    24  
     24
     25
     26
     27
     28
    2529  <!-- Define domains and groups of domains -->
    2630  <domain_definition>
    27     <domain id="dom_glo" data_dim="2" />
     31
     32    <domain id="dom_glo"  data_dim="1" />
     33
     34    <domain id="greordered"  domain_ref="dom_glo">
     35      <reorder_domain invert_lat="true" shift_lon_fraction="0.5" min_lon="0" max_lon="360" />
     36    </domain>
     37   
     38    <domain id="dom_regular" ni_glo="144" nj_glo="142" type="rectilinear"  >
     39      <generate_rectilinear_domain/>
     40      <interpolate_domain order="1"/>
     41    </domain>
     42   
     43    <domain id="dom_out" domain_ref="dom_glo"/>
     44   
    2845  </domain_definition>
    2946 
    … …  
    3653  <grid_definition>
    3754    <grid id="grid_scalar" >
    38     <scalar/>
     55      <scalar/>
    3956    </grid>
    4057  </grid_definition>
    … …  
    4259  <!-- Define groups of vertical axes -->
    4360  <axis_definition>
    44     <axis id="presnivs" standard_name="Vertical levels" unit="Pa">
     61    <axis id="time_month" n_glo="12" value="(0,11) [1 2 3 4 5 6 7 8 9 10 11 12]"/>
     62    <axis id="time_year" unit="day" />
     63    <axis id="presnivs" standard_name="Vertical levels" unit="Pa"/>
     64    <axis id="Ahyb" standard_name="Ahyb comp of Hyb Cord" unit="Pa"/>
     65    <axis id="Bhyb" standard_name="Bhyb comp of Hyb Cord" unit=""/>
     66    <axis id="Ahyb_inter" standard_name="A comp of Hyb Cord at interface" unit="Pa"/>
     67    <axis id="Bhyb_inter" standard_name="B comp of Hyb Cord at interface" unit=""/>
     68    <axis id="Alt" standard_name="Height approx for scale heigh of 8km at levels" unit="km"/>
     69    <axis id="plev" standard_name="model_level_number" unit="Pa"/>
     70    <axis id="klev"  prec="8" long_name="number of layers" standard_name="number of layers" unit="1" />
     71    <axis id="klevp1"  prec="8" long_name="number of layer interfaces" standard_name="number of layer interfaces" unit="1" />
     72    <axis id="bnds" standard_name="bounds" unit="1" />
     73    <axis id="spectband" standard_name="Sensor Band Central Radiation Wavenumber" unit="m-1"/>
     74    <axis id="axis_lat" standard_name="Latitude axis">
     75      <reduce_domain operation="average" direction="iDir" />
    4576    </axis>
    46     <axis id="Ahyb" standard_name="Ahyb comp of Hyb Cord" unit="Pa">
    47     </axis>
    48     <axis id="Bhyb" standard_name="Bhyb comp of Hyb Cord" unit="">
    49     </axis>
    50     <axis id="Ahyb_inter" standard_name="A comp of Hyb Cord at interface" unit="Pa">
    51     </axis>
    52     <axis id="Bhyb_inter" standard_name="B comp of Hyb Cord at interface" unit="">
    53     </axis>
    54     <axis id="Alt" standard_name="Height approx for scale heigh of 8km at levels" unit="km">
    55     </axis>
    56     <axis id="plev" standard_name="model_level_number" unit="Pa">
    57     </axis>
    58     <axis id="klev"  prec="8" long_name = "hybrid sigma pressure coordinate"
    59           standard_name ="atmosphere_hybrid_sigma_pressure_coordinate" unit="1">
    60     </axis>
    61     <axis id="bnds" standard_name="bounds" unit="1" >
    62     </axis>
    6377
    6478<!-- Cosp axis definitions-->
    65     <axis id="height" standard_name="Cosp levels" unit="m">
    66     </axis>
    67     <axis id="height_mlev" standard_name="height_mlev" unit="m">
    68     </axis>
    69     <axis id="sza" standard_name="solar_zenith_angle" unit="degrees">
    70     </axis>
    71     <axis id="pressure2" standard_name="pressure" unit="mb">
    72     </axis>
    73     <axis id="column" standard_name="column" unit="count">
    74     </axis>
    75     <axis id="temp" standard_name="temperature" unit="K">
    76     </axis>
    77     <axis id="cth16" name="cth" standard_name="altitude" unit="m">
    78     </axis>
    79     <axis id="ReffIce" standard_name="ReffIce" unit="microne" >
    80     </axis>
    81     <axis id="ReffLiq" standard_name="ReffLiq" unit="microne" >
    82     </axis>
    83     <axis id="scatratio" standard_name="scatratio" unit="1" >
    84     </axis>
    85     <axis id="dbze" standard_name="dbze" unit="dBZ" >
    86     </axis>
    87     <axis id="tau" standard_name="tau" unit="1" >
    88     </axis>
     79    <axis id="height" standard_name="Cosp levels" unit="m"/>
     80    <axis id="height_mlev" standard_name="height_mlev" unit="m"/>
     81    <axis id="sza" standard_name="solar_zenith_angle" unit="degrees"/>
     82    <axis id="pressure2" standard_name="pressure" unit="mb"/>
     83    <axis id="column" standard_name="column" unit="count"/>
     84    <axis id="temp" standard_name="temperature" unit="K"/>
     85    <axis id="cth16" standard_name="altitude" unit="m"/>
     86    <axis id="ReffIce" standard_name="ReffIce" unit="microne" />
     87    <axis id="ReffLiq" standard_name="ReffLiq" unit="microne" />
     88    <axis id="scatratio" standard_name="scatratio" unit="1" />
     89    <axis id="dbze" standard_name="dbze" unit="dBZ" />
     90    <axis id="tau" standard_name="tau" unit="1" />             
    8991  </axis_definition>
    9092
    9193  <grid_definition>
    92 <!-- Define Scalar grid for GHG, orbital parameters and solar constants -->
    93     <grid id="grid_scalar">
    94     </grid>
    9594
    9695    <grid id="klev_bnds"> <axis axis_ref="klev" /> <axis axis_ref="bnds" /> </grid>
     96    <grid id="klevp1_bnds"> <axis axis_ref="klevp1" /> <axis axis_ref="bnds" /> </grid>
     97
    9798
    9899     <grid id="grid_glo">
    99         <domain id="dom_glo" />
    100      </grid>
     100        <domain domain_ref="dom_glo" />
     101     </grid>
     102
     103     <grid id="grid_out">
     104        <domain domain_ref="dom_out" />
     105     </grid>
     106
    101107
    102108     <grid id="grid_glo_presnivs">
    103         <domain id="dom_glo" />
    104         <axis id="presnivs" />
     109        <domain domain_ref="dom_glo" />
     110        <axis axis_ref="presnivs" />
     111     </grid>
     112
     113     <grid id="grid_out_presnivs">
     114        <domain domain_ref="dom_out" />
     115        <axis axis_ref="presnivs" />
    105116     </grid>
    106117
    107118
    108119     <grid id="grid_glo_plev">
    109         <domain id="dom_glo" />
    110         <axis id="plev" />
    111      </grid>
     120        <domain domain_ref="dom_glo" />
     121        <axis axis_ref="plev" />
     122     </grid>
     123
     124     <grid id="grid_out_plev">
     125        <domain domain_ref="dom_out" />
     126        <axis axis_ref="plev" />
     127     </grid>
     128
     129
     130     <grid id="grid_glo_spectband">
     131        <domain domain_ref="dom_glo" />
     132        <axis axis_ref="spectband" />
     133     </grid>
     134
     135     <grid id="grid_out_spectband">
     136        <domain domain_ref="dom_out" />
     137        <axis axis_ref="spectband" />
     138     </grid>
     139     
    112140
    113141     <grid id="grid_glo_height">
    114         <domain id="dom_glo" />
    115         <axis id="height" />
    116      </grid>
     142        <domain domain_ref="dom_glo" />
     143        <axis axis_ref="height" />
     144     </grid>
     145     <grid id="grid_out_height">
     146        <domain domain_ref="dom_out" />
     147        <axis axis_ref="height" />
     148     </grid>
     149
    117150
    118151     <grid id="grid_glo_heightmlev">
    119         <domain id="dom_glo" />
    120         <axis id="height_mlev" />
    121      </grid>
     152        <domain domain_ref="dom_glo" />
     153        <axis axis_ref="height_mlev" />
     154     </grid>
     155     <grid id="grid_out_heightmlev">
     156        <domain domain_ref="dom_out" />
     157        <axis axis_ref="height_mlev" />
     158     </grid>
     159
    122160
    123161     <grid id="grid_glo_temp">
    124         <domain id="dom_glo" />
    125         <axis id="temp" />
    126      </grid>
     162        <domain domain_ref="dom_glo" />
     163        <axis axis_ref="temp" />
     164     </grid>
     165     <grid id="grid_out_temp">
     166        <domain domain_ref="dom_out" />
     167        <axis axis_ref="temp" />
     168     </grid>
     169
     170
    127171
    128172     <grid id="grid_glo_sza">
    129         <domain id="dom_glo" />
    130         <axis id="sza" />
    131      </grid>
     173        <domain domain_ref="dom_glo" />
     174        <axis axis_ref="sza" />
     175     </grid>
     176     <grid id="grid_out_sza">
     177        <domain domain_ref="dom_out" />
     178        <axis axis_ref="sza" />
     179     </grid>
     180
    132181
    133182     <grid id="grid_glo_column">
    134         <domain id="dom_glo" />
    135         <axis id="column" />
     183        <domain domain_ref="dom_glo" />
     184        <axis axis_ref="column" />
     185     </grid>
     186     <grid id="grid_out_column">
     187        <domain domain_ref="dom_out" />
     188        <axis axis_ref="column" />
    136189     </grid>
    137190
    … …  
    139192<!-- Define 4D grids for Cosp simulator -->
    140193     <grid id="grid4Dcol">
    141         <domain id="dom_glo" />
    142         <axis id="height_mlev" />
    143         <axis id="column" />
     194        <domain domain_ref="dom_glo" />
     195        <axis axis_ref="column" />
     196        <axis axis_ref="height_mlev" />
     197     </grid>
     198     <grid id="grid4Dcol_out">
     199        <domain domain_ref="dom_out" />
     200        <axis axis_ref="column" />
     201        <axis axis_ref="height_mlev" />
    144202     </grid>
    145203
    146204     <grid id="grid4Dsrbin">
    147         <domain id="dom_glo" />
    148         <axis id="height" />
    149         <axis id="scatratio" />
    150      </grid>
     205        <domain domain_ref="dom_glo" />
     206        <axis axis_ref="height" />
     207        <axis axis_ref="scatratio" />
     208     </grid>
     209     <grid id="grid4Dsrbin_out">
     210        <domain domain_ref="dom_out" />
     211        <axis axis_ref="height" />
     212        <axis axis_ref="scatratio" />
     213     </grid>
     214
    151215
    152216     <grid id="grid4Ddbze">
    153         <domain id="dom_glo" />
    154         <axis id="height" />
    155         <axis id="dbze" />
     217        <domain domain_ref="dom_glo" />
     218        <axis axis_ref="height" />
     219        <axis axis_ref="dbze" />
     220     </grid>
     221     <grid id="grid4Ddbze_out">
     222        <domain domain_ref="dom_out" />
     223        <axis axis_ref="height" />
     224        <axis axis_ref="dbze" />
    156225     </grid>
    157226
    158227     <grid id="grid4Dtau">
    159         <domain id="dom_glo" />
    160         <axis id="tau" />
    161         <axis id="pressure2" />
     228        <domain domain_ref="dom_glo" />
     229        <axis axis_ref="tau" />
     230        <axis axis_ref="pressure2" />
     231     </grid>
     232     <grid id="grid4Dtau_out">
     233        <domain domain_ref="dom_out" />
     234        <axis axis_ref="tau" />
     235        <axis axis_ref="pressure2" />
    162236     </grid>
    163237
    164238     <grid id="grid4Dmisr">
    165         <domain id="dom_glo" />
    166         <axis id="tau" />
    167         <axis id="cth16" />
     239        <domain domain_ref="dom_glo" />
     240        <axis axis_ref="cth16" />
     241        <axis axis_ref="tau" />
     242     </grid>
     243     <grid id="grid4Dmisr_out">
     244        <domain domain_ref="dom_out" />
     245        <axis axis_ref="cth16" />
     246        <axis axis_ref="tau" />
    168247     </grid>
    169248
    170249     <grid id="grid4Dreffi">
    171         <domain id="dom_glo" />
    172         <axis id="tau" />
    173         <axis id="ReffIce" />
    174      </grid>
    175 
     250        <domain domain_ref="dom_glo" />
     251        <axis axis_ref="tau" />
     252        <axis axis_ref="ReffIce" />
     253     </grid>
     254     <grid id="grid4Dreffi_out">
     255        <domain domain_ref="dom_out" />
     256        <axis axis_ref="tau" />
     257        <axis axis_ref="ReffIce" />
     258     </grid>
     259
     260     <grid id="grid4Dreffl_out">
     261        <domain domain_ref="dom_glo" />
     262        <axis axis_ref="tau" />
     263        <axis axis_ref="ReffLiq" />
     264     </grid>
    176265     <grid id="grid4Dreffl">
    177         <domain id="dom_glo" />
    178         <axis id="tau" />
    179         <axis id="ReffLiq" />
    180      </grid>
     266        <domain domain_ref="dom_out" />
     267        <axis axis_ref="tau" />
     268        <axis axis_ref="ReffLiq" />
     269     </grid>
     270
     271
    181272     <grid id="grid4Dcol2">
    182         <domain id="dom_glo" />
    183         <axis id="height" />
    184         <axis id="column" />
    185      </grid>
     273        <domain domain_ref="dom_glo" />
     274        <axis axis_ref="height" />
     275        <axis axis_ref="column" />
     276     </grid>
     277     <grid id="grid4Dcol2_out">
     278        <domain domain_ref="dom_out" />
     279        <axis axis_ref="height" />
     280        <axis axis_ref="column" />
     281     </grid>
     282
     283 <!-- Grid definitions to allow summing of a 3D variable -->   
     284      <grid id="grid_3D" >
     285         <domain domain_ref="dom_glo" />
     286         <axis axis_ref="lev" />
     287       </grid>
     288
     289       <grid id="grid_sum_axis">
     290         <domain domain_ref="dom_glo" />
     291         <scalar>
     292            <reduce_axis operation="sum" />
     293         </scalar>
     294       </grid>
     295
     296       <grid id="grid_sum">
     297         <scalar >
     298           <reduce_domain operation="sum" local="true" />
     299           <reduce_scalar operation="sum" />
     300         </scalar>
     301         <scalar/>
     302       </grid>
    186303
    187304  </grid_definition>

  • LMDZ6/branches/DYNAMICO-conv/bld.cfg

    r3411 r3413  
    6666# Example of how to set specific compiling options for a specific file
    6767# -> this can be including in the arch.opt file
    68 #bld::tool::fflags::phys::readaerosol         %BASE_FFLAGS %PROD_FFLAGS  %INCDIR -C hopt -pi auto
     68#bld::tool::fflags::phys::physiq_mod         %BASE_FFLAGS %DEBUG_FFLAGS  %INCDIR
    6969
    7070

  • LMDZ6/branches/DYNAMICO-conv/build_gcm

    r1907 r3413  
    4545# run "fcm build" command
    4646${dirname}fcm build $*
     47status=$?
    4748
    4849# cleanup
    4950\rm -f '.lock'
    5051
     52if [ $status -ne 0 ]; then
     53  exit $status
     54fi
     55

  • LMDZ6/branches/DYNAMICO-conv/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90

    r2941 r3413  
    100100
    101101  deg2rad = pi/180.0
    102 
     102  y(:,:,:)=0  !ym warning unitialized variable
     103 
    103104! Compute psol AND tsol, knowing phis.
    104105!*******************************************************************************

  • LMDZ6/branches/DYNAMICO-conv/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90

    r2899 r3413  
    4444    zmax0,fevap, rnebcon,falb_dir, wake_fip,    agesno,  detr_therm, pbl_tke,  &
    4545    phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
    46     prw_ancien
     46    prw_ancien, sollwdown
    4747  USE comconst_mod, ONLY: pi, dtvr
    4848
    … …  
    200200  solsw      = 165.
    201201  sollw      = -53.
     202!ym warning missing init for sollwdown => set to 0
     203  sollwdown  = 0.
    202204  t_ancien   = 273.15
    203205  q_ancien   = 0.
    … …  
    314316  ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation
    315317  ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy
     318  zsig0(:,:)=0   !ym uninitialized variable
     319  zgam0(:,:)=0   !ym uninitialized variable
    316320  ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
     321  zthe0(:,:)=0   !ym uninitialized variable
    317322  ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
    318323

  • LMDZ6/branches/DYNAMICO-conv/libf/misc/wxios.F90

    r3411 r3413  
    152152!$OMP MASTER
    153153        !Initialisation du contexte:
     154        !!CALL xios_context_initialize(g_ctx_name, g_comm)
    154155        CALL xios_context_initialize(g_ctx_name, COMM_LMDZ_PHY)
    155156        CALL xios_get_handle(g_ctx_name, xios_ctx)    !Récupération
    … …  
    301302        CALL xios_set_domain_attr_hdl(dom, nj_glo=nbp_lat, jbegin=jj_begin-1, nj=jj_nb, data_dim=2)
    302303        CALL xios_set_domain_attr_hdl(dom, lonvalue_1d=io_lon(1:nbp_lon), latvalue_1d=io_lat(jj_begin:jj_end))
     304        CALL xios_set_domain_attr("dom_out", domain_ref=dom_id)
    303305
    304306        IF (.NOT.is_sequential) THEN
    … …  
    332334        USE mod_phys_lmdz_para
    333335        USE nrtype, ONLY : PI
     336        USE ioipsl_getin_p_mod, ONLY : getin_p
    334337        IMPLICIT NONE
    335338        CHARACTER(len=*),INTENT(IN) :: dom_id ! domain identifier
    … …  
    339342        REAL :: boundslat_mpi(klon_mpi,nvertex)
    340343        INTEGER :: ind_cell_glo_mpi(klon_mpi)
    341         TYPE(xios_domaingroup) :: dom
    342 
     344        TYPE(xios_domain) :: dom
     345        LOGICAL :: remap_output
    343346
    344347        CALL gather_omp(longitude*180/PI,lon_mpi)
    … …  
    348351        CALL gather_omp(ind_cell_glo,ind_cell_glo_mpi)
    349352       
     353        remap_output=.TRUE.
     354        CALL getin_p("remap_output",remap_output)
    350355
    351356!$OMP MASTER
    352         CALL xios_get_domaingroup_handle(dom_id, dom)
     357        CALL xios_get_domain_handle(dom_id, dom)
    353358       
    354359        !On parametrise le domaine:
    … …  
    357362                           bounds_lon_1d=TRANSPOSE(boundslon_mpi), bounds_lat_1d=TRANSPOSE(boundslat_mpi) )
    358363        CALL xios_set_attr(dom, i_index=ind_cell_glo_mpi(:)-1)
     364        IF (remap_output) THEN
     365          CALL xios_set_domain_attr("dom_out", domain_ref="dom_regular")
     366          CALL xios_set_fieldgroup_attr("dom_out", domain_ref="dom_regular")
     367          CALL xios_set_fieldgroup_attr("remap_expr", expr="@this_ref")
     368        ENDIF
    359369!$OMP END MASTER
    360370

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/conf_phys_m.F90

    r3411 r3413  
    10091009    ! - 1 = stratospheric aerosols scaled from 550 nm AOD
    10101010    ! - 2 = stratospheric aerosol properties from CMIP6
     1011    !Option 2 is only available with RRTM, this is tested later on
    10111012    !Config Def  = 0
    10121013    !Config Help = Used in physiq.F
    … …  
    10141015    !
    10151016    flag_aerosol_strat_omp = 0
    1016     IF (iflag_rrtm_omp==1) THEN
    1017       CALL getin('flag_aerosol_strat',flag_aerosol_strat_omp)
    1018     ENDIF
     1017    CALL getin('flag_aerosol_strat',flag_aerosol_strat_omp)
    10191018
    10201019    !

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/cosp/phys_cosp.F90

    r3411 r3413  
    330330      !$OMP END MASTER
    331331      !$OMP BARRIER
    332         debut_cosp=.false.
    333332      endif ! debut_cosp
    334333!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    … …  
    341340!        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
    342341!#else
    343        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
     342     if (.NOT. debut_cosp) call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
    344343!#endif
    345344!!
    … …  
    350349
    351350!       print *, 'Calling write output'
    352        call cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, missing_val, &
     351     if (.NOT. debut_cosp) call cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, missing_val, &
    353352                               cfg, gbx, vgrid, sglidar, sgradar, stlidar, stradar, &
    354353                               isccp, misr, modis)
    … …  
    376375!  call system_clock(t1,count_rate,count_max)
    377376!  print *,(t1-t0)*1.0/count_rate
     377    if (debut_cosp) then
     378      debut_cosp=.false.
     379    endif
    378380 
    379381  CONTAINS

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_unstruct.F90

    r3323 r3413  
    135135    rain_fall = 0.; snow_fall = 0.
    136136    solsw = 165.;   sollw = -53.
     137!ym warning missing init for sollwdown => set to 0
     138  sollwdown  = 0.
     139   
     140   
    137141    t_ancien = 273.15
    138142    u_ancien=0

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/cv3_routines.F90

    r3411 r3413  
    37563756        IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
    37573757          qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
    3758 !jyg<   Bug correction 20180620
    3759 !      PROBLEM: Should not qent(il,i,i) be taken into account even if nent(il,i)/=0?
    3760 !!        qtment(il, i) = qent(il,k,i) + qtment(il,i)
    3761 ! cld
    3762           qtment(il, i) = qent(il,i,i) + qtment(il,i)
    3763 ! cld
    3764 !>jyg
    3765 
     3758          qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
    37663759          nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
    37673760        END IF                                                                   ! cld

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/limit_read_mod.F90

    r3336 r3413  
    258258          ENDIF
    259259          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
    260           WRITE(abort_message,'(a,2(i3,a))')'limit.nc records number (',nn,') does no'//&
     260          WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
    261261            't match year length (',year_len,')'
    262262          IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
    … …  
    336336
    337337!$OMP MASTER  ! Only master thread
    338        IF (is_mpi_root) THEN ! Only master processus
     338       IF (is_mpi_root) THEN ! Only master processus!
    339339
    340340          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_cal_mod.F90

    r3336 r3413  
    3636 
    3737  SUBROUTINE phys_cal_init(annee_ref,day_ref)
    38   USE IOIPSL, ONLY:  ymds2ju, ioconf_calendar
    39   USE mod_phys_lmdz_para, ONLY:  is_master,is_omp_master
    40   USE ioipsl_getin_p_mod, ONLY: getin_p
    41   IMPLICIT NONE
     38
     39    USE IOIPSL, ONLY:  ymds2ju, ioconf_calendar
     40    USE mod_phys_lmdz_para, ONLY:  is_master,is_omp_master
     41    USE ioipsl_getin_p_mod, ONLY: getin_p
     42
     43    IMPLICIT NONE
    4244    INTEGER,INTENT(IN) :: annee_ref
    4345    INTEGER,INTENT(IN) :: day_ref
    … …  
    4648    calend = 'earth_360d' ! default
    4749    CALL getin_p("calend",calend)
    48      
     50
    4951    IF (is_omp_master) THEN
    5052      IF (calend == 'earth_360d') THEN
    … …  
    5961    ENDIF
    6062!$OMP BARRIER
    61 
     63     
    6264    CALL ymds2ju(annee_ref, 1, day_ref, 0., jD_ref)
    6365    jD_ref=INT(jD_ref)

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_state_var_mod.F90

    r3411 r3413  
    2323!$OMP THREADPRIVATE(cvpas)
    2424!$OMP THREADPRIVATE(wkpas)
    25       REAL, SAVE :: phys_tstep, solaire_etat0
     25      REAL, SAVE :: phys_tstep=0, solaire_etat0
    2626!$OMP THREADPRIVATE(phys_tstep, solaire_etat0)
    2727
    … …  
    607607        du_gwd_rando(:,:)=0.
    608608      ENDIF
    609       if (.not. ok_hines .and. ok_gwd_rando) allocate(du_gwd_front(klon, klev))
    610 
     609      IF (.not. ok_hines .and. ok_gwd_rando) THEN
     610        ALLOCATE(du_gwd_front(klon, klev))
     611        du_gwd_front(:,:) = 0 !ym missing init   
     612      ENDIF
    611613END SUBROUTINE phys_state_var_init
    612614

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90

    r3411 r3413  
    258258    USE climoz_mod
    259259    USE limit_read_mod, ONLY : init_limit_read
     260    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
     261    USE readaerosol_mod, ONLY : init_aero_fromfile
     262    USE readaerosolstrato_m, ONLY : init_readaerosolstrato
    260263
    261264
    … …  
    11951198       CALL phys_state_var_init(read_climoz)
    11961199       CALL phys_output_var_init
     1200       IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
     1201
    11971202       print*, '================================================='
    11981203       !
    … …  
    14921497       !
    14931498
    1494        CALL create_climoz(read_climoz)
     1499!       CALL create_climoz(read_climoz)
     1500       CALL init_aero_fromfile(flag_aerosol)  !! initialise aero from file for XIOS interpolation (unstructured_grid)
     1501       CALL init_readaerosolstrato(flag_aerosol_strat)  !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
     1502
     1503#ifdef CPP_COSP
     1504       CALL phys_cosp(itap,phys_tstep,freq_cosp, &
     1505               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
     1506               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
     1507               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
     1508               JrNt,ref_liq,ref_ice, &
     1509               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
     1510               zu10m,zv10m,pphis, &
     1511               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
     1512               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
     1513               prfl(:,1:klev),psfl(:,1:klev), &
     1514               pmflxr(:,1:klev),pmflxs(:,1:klev), &
     1515               mr_ozone,cldtau, cldemi)
     1516#endif
    14951517
    14961518       CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
    … …  
    15041526       IF (is_omp_master) CALL xios_update_calendar(1)
    15051527#endif
     1528       IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
    15061529       CALL create_etat0_limit_unstruct
    1507 
    15081530       CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
    15091531
    … …  
    15351557       ENDIF
    15361558
    1537        CALL printflag( tabcntr0,radpas,ok_journe, &
    1538             ok_instan, ok_region )
    1539        !
    15401559!       IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
    15411560!          WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
    … …  
    47334752       !         write(97) u_seri,v_seri,t_seri,q_seri
    47344753       !         close(97)
    4735 !       !$OMP MASTER
    4736        IF (read_climoz >= 1) THEN
    4737           IF (is_mpi_root) THEN
    4738              CALL nf95_close(ncid_climoz)
    4739           ENDIF
    4740           DEALLOCATE(press_edg_climoz) ! pointer
    4741           DEALLOCATE(press_cen_climoz) ! pointer
    4742        ENDIF
    4743 !       !$OMP END MASTER
     4754     
     4755       IF (is_omp_master) THEN
     4756       
     4757         IF (read_climoz >= 1) THEN
     4758           IF (is_mpi_root) CALL nf95_close(ncid_climoz)
     4759            DEALLOCATE(press_edg_climoz) ! pointer
     4760            DEALLOCATE(press_cen_climoz) ! pointer
     4761         ENDIF
     4762       
     4763       ENDIF
    47444764#ifdef CPP_XIOS
    47454765       IF (is_omp_master) CALL xios_context_finalize

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/radlwsw_m.F90

    r3411 r3413  
    385385  lldebug=.FALSE.
    386386
     387  ztopsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
     388  ztopsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
     389  zsolsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
     390  zsolsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    387391 
    388392  !
    … …  
    597601      ztopswaiaero(i)=0.
    598602      zsolswaiaero(i)=0.
    599       ztopsw_aero(i,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    600       ztopsw0_aero(i,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    601       zsolsw_aero(i,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    602       zsolsw0_aero(i,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
    603603      ENDDO
    604604!     print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosol.F90

    r2841 r3413  
    44
    55  REAL, SAVE :: not_valid=-333.
     6 
     7  INTEGER, SAVE :: nbp_lon_src
     8!$OMP THREADPRIVATE(nbp_lon_src) 
     9  INTEGER, SAVE :: nbp_lat_src
     10!$OMP THREADPRIVATE(nbp_lat_src) 
     11  REAL, ALLOCATABLE, SAVE    :: psurf_interp(:,:)
     12!$OMP THREADPRIVATE(psurf_interp) 
    613
    714CONTAINS
    … …  
    167174
    168175
     176  SUBROUTINE init_aero_fromfile(flag_aerosol)
     177  USE netcdf
     178  USE mod_phys_lmdz_para
     179  USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
     180  USE xios
     181  IMPLICIT NONE
     182  INTEGER, INTENT(IN) :: flag_aerosol
     183  REAL,ALLOCATABLE :: lat_src(:)
     184  REAL,ALLOCATABLE :: lon_src(:)
     185  CHARACTER(LEN=*),PARAMETER :: file_aerosol='aerosols.nat.nc'
     186  CHARACTER(LEN=*),PARAMETER :: file_so4='so4.nat.nc'
     187  INTEGER :: klev_src
     188  INTEGER :: ierr ,ncid, dimID, varid
     189  REAL :: null_array(0)
     190
     191  IF (flag_aerosol>0 .AND. grid_type==unstructured) THEN
     192 
     193    IF (is_omp_root) THEN
     194 
     195      IF (is_mpi_root) THEN
     196   
     197        IF (nf90_open(TRIM(file_aerosol), NF90_NOWRITE, ncid) /= NF90_NOERR) THEN
     198          CALL check_err( nf90_open(TRIM(file_so4), NF90_NOWRITE, ncid), "pb open "//trim(file_so4) )
     199        ENDIF
     200
     201        ! Read and test longitudes
     202        CALL check_err( nf90_inq_dimid(ncid, "lon", dimID),"pb inq dim lon")
     203        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src),"pb inq dim lon")
     204        CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
     205        ALLOCATE(lon_src(nbp_lon_src))
     206        CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
     207
     208        ! Read and test latitudes
     209        CALL check_err( nf90_inq_dimid(ncid, "lat", dimID),"pb inq dim lat")
     210        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src),"pb inq dim lat")
     211        CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
     212        ALLOCATE(lat_src(nbp_lat_src))
     213        CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
     214        IF (nf90_inq_dimid(ncid, 'lev', dimid) /= NF90_NOERR) THEN
     215          IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= NF90_NOERR) THEN
     216             CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid),'dimension lev,PRESNIVS or presnivs not in file')
     217          ENDIF
     218        ENDIF
     219        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
     220        CALL check_err( nf90_close(ncid),"pb in close" )   
     221      ENDIF
     222
     223      CALL bcast_mpi(nbp_lat_src)
     224      CALL bcast_mpi(nbp_lon_src)
     225      CALL bcast_mpi(klev_src)
     226
     227      IF (is_mpi_root ) THEN
     228        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=nbp_lat_src, jbegin=0, latvalue_1d=lat_src)
     229        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=nbp_lon_src, ibegin=0, lonvalue_1d=lon_src)
     230      ELSE
     231        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=0, jbegin=0, latvalue_1d=null_array )
     232        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=0, ibegin=0, lonvalue_1d=null_array)
     233      ENDIF
     234      CALL xios_set_axis_attr("axis_aerosol",n_glo=klev_src)
     235      CALL xios_set_fieldgroup_attr("aerosols", enabled=.TRUE.)
     236 
     237    ENDIF
     238   
     239  ENDIF   
     240 
     241  END SUBROUTINE init_aero_fromfile
     242
     243
     244   
     245
    169246  SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
    170247!****************************************************************************************
    … …  
    187264    USE netcdf
    188265    USE dimphy
    189     USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, &
    190                                  grid2Dto1D_glo
     266    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, &
     267                                 grid2Dto1D_glo, grid_type, unstructured
    191268    USE mod_phys_lmdz_para
    192269    USE iophy, ONLY : io_lon, io_lat
    193270    USE print_control_mod, ONLY: lunout
     271    USE xios
    194272
    195273    IMPLICIT NONE
    … …  
    205283    REAL, POINTER, DIMENSION(:)           :: pt_b         ! Pointer for describing the vertical levels     
    206284    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
     285    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year_mpi  ! Pointer-variabale from file, 12 month, grid : klon,klev_src
    207286    REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
     287    REAL, DIMENSION(klon_mpi,12)          :: psurf_out_mpi    ! Surface pression for 12 months
    208288    REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
     289    REAL, DIMENSION(klon_mpi,12)          :: load_out_mpi     ! Aerosol mass load in each column
    209290    INTEGER                               :: nbr_tsteps   ! number of month in file read
    210291
    … …  
    220301    REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
    221302
    222     REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
     303    REAL,  ALLOCATABLE                    :: psurf_glo2D(:,:,:)   ! Surface pression for 12 months on dynamics global grid
    223304    REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
    224     REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: load_glo2D    ! Load for 12 months on dynamics global grid
     305    REAL,  ALLOCATABLE                    :: load_glo2D(:,:,:)    ! Load for 12 months on dynamics global grid
    225306    REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
    226     REAL, DIMENSION(nbp_lon,nbp_lat)      :: vartmp
    227     REAL, DIMENSION(nbp_lon)              :: lon_src              ! longitudes in file
    228     REAL, DIMENSION(nbp_lat)              :: lat_src, lat_src_inv ! latitudes in file
     307    REAL, ALLOCATABLE, DIMENSION(:,:)     :: vartmp
     308    REAL, ALLOCATABLE,DIMENSION(:)        :: lon_src              ! longitudes in file
     309    REAL, ALLOCATABLE,DIMENSION(:)        :: lat_src, lat_src_inv ! latitudes in file
    229310    LOGICAL                               :: new_file             ! true if new file format detected
    230311    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
    231 
    232 
     312    INTEGER                               :: nbp_lon, nbp_lat
     313    LOGICAL,SAVE                          :: first=.TRUE.
     314!$OMP THREADPRIVATE(first)
     315   
     316    IF (grid_type==unstructured) THEN
     317      nbp_lon=nbp_lon_src
     318      nbp_lat=nbp_lat_src
     319    ELSE
     320      nbp_lon=nbp_lon_
     321      nbp_lat=nbp_lat_
     322    ENDIF
     323   
     324    IF (is_mpi_root) THEN
     325   
     326      ALLOCATE(psurf_glo2D(nbp_lon,nbp_lat,12))
     327      ALLOCATE(load_glo2D(nbp_lon,nbp_lat,12))
     328      ALLOCATE(vartmp(nbp_lon,nbp_lat))
     329      ALLOCATE(lon_src(nbp_lon))
     330      ALLOCATE(lat_src(nbp_lat))
     331      ALLOCATE(lat_src_inv(nbp_lat))
     332    ELSE
     333      ALLOCATE(varyear(0,0,0,0))
     334      ALLOCATE(psurf_glo2D(0,0,0))
     335      ALLOCATE(load_glo2D(0,0,0))
     336    ENDIF
     337           
    233338    ! Deallocate pointers
    234339    IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
    … …  
    245350       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) )
    246351
     352
     353       IF (grid_type/=unstructured) THEN
     354
    247355! Test for equal longitudes and latitudes in file and model
    248356!****************************************************************************************
    249357       ! Read and test longitudes
    250        CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
    251        CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
     358         CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
     359         CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
    252360       
    253        IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
    254           WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
    255           WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
    256           WRITE(lunout,*) 'longitudes in model :', io_lon
     361         IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
     362            WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
     363            WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
     364            WRITE(lunout,*) 'longitudes in model :', io_lon
    257365         
    258           CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
    259        END IF
    260 
    261        ! Read and test latitudes
    262        CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
    263        CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
    264 
    265        ! Invert source latitudes
    266        DO j = 1, nbp_lat
    267           lat_src_inv(j) = lat_src(nbp_lat +1 -j)
    268        END DO
    269 
    270        IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
    271           ! Latitudes are the same
    272           invert_lat=.FALSE.
    273        ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
    274           ! Inverted source latitudes correspond to model latitudes
    275           WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
    276           invert_lat=.TRUE.
    277        ELSE
    278           WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
    279           WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
    280           WRITE(lunout,*) 'latitudes in model :', io_lat
    281           CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
    282        END IF
    283 
     366            CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
     367         END IF
     368
     369         ! Read and test latitudes
     370         CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
     371         CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
     372
     373         ! Invert source latitudes
     374         DO j = 1, nbp_lat
     375            lat_src_inv(j) = lat_src(nbp_lat +1 -j)
     376         END DO
     377
     378         IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
     379            ! Latitudes are the same
     380            invert_lat=.FALSE.
     381         ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
     382            ! Inverted source latitudes correspond to model latitudes
     383            WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
     384            invert_lat=.TRUE.
     385         ELSE
     386            WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
     387            WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
     388            WRITE(lunout,*) 'latitudes in model :', io_lat
     389            CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
     390         END IF
     391       ENDIF
    284392
    285393! 2) Check if old or new file is avalabale.
    … …  
    487595       END IF
    488596
    489 !     - Invert latitudes if necessary
    490        DO imth=1, 12
    491           IF (invert_lat) THEN
    492 
    493              ! Invert latitudes for the variable
    494              varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
    495              DO k=1,klev_src
    496                 DO j=1,nbp_lat
    497                    DO i=1,nbp_lon
    498                       varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
    499                    END DO
    500                 END DO
    501              END DO
     597
     598       IF (grid_type/=unstructured) THEN
     599  !     - Invert latitudes if necessary
     600         DO imth=1, 12
     601            IF (invert_lat) THEN
     602
     603               ! Invert latitudes for the variable
     604               varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
     605               DO k=1,klev_src
     606                  DO j=1,nbp_lat
     607                     DO i=1,nbp_lon
     608                        varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
     609                     END DO
     610                  END DO
     611               END DO
    502612             
    503              ! Invert latitudes for surface pressure
    504              vartmp(:,:) = psurf_glo2D(:,:,imth)
    505              DO j=1,nbp_lat
    506                 DO i=1,nbp_lon
    507                    psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
    508                 END DO
    509              END DO
     613               ! Invert latitudes for surface pressure
     614               vartmp(:,:) = psurf_glo2D(:,:,imth)
     615               DO j=1,nbp_lat
     616                  DO i=1,nbp_lon
     617                     psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
     618                  END DO
     619               END DO
    510620             
    511              ! Invert latitudes for the load
    512              vartmp(:,:) = load_glo2D(:,:,imth)
    513              DO j=1,nbp_lat
    514                 DO i=1,nbp_lon
    515                    load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
    516                 END DO
    517              END DO
    518           END IF ! invert_lat
     621               ! Invert latitudes for the load
     622               vartmp(:,:) = load_glo2D(:,:,imth)
     623               DO j=1,nbp_lat
     624                  DO i=1,nbp_lon
     625                     load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
     626                  END DO
     627               END DO
     628            END IF ! invert_lat
    519629             
    520           ! Do zonal mead at poles and distribut at whole first and last latitude
    521           DO k=1, klev_src
    522              npole=0.  ! North pole, j=1
    523              spole=0.  ! South pole, j=nbp_lat       
    524              DO i=1,nbp_lon
    525                 npole = npole + varyear(i,1,k,imth)
    526                 spole = spole + varyear(i,nbp_lat,k,imth)
    527              END DO
    528              npole = npole/REAL(nbp_lon)
    529              spole = spole/REAL(nbp_lon)
    530              varyear(:,1,    k,imth) = npole
    531              varyear(:,nbp_lat,k,imth) = spole
    532           END DO
    533        END DO ! imth
     630            ! Do zonal mead at poles and distribut at whole first and last latitude
     631            DO k=1, klev_src
     632               npole=0.  ! North pole, j=1
     633               spole=0.  ! South pole, j=nbp_lat       
     634               DO i=1,nbp_lon
     635                  npole = npole + varyear(i,1,k,imth)
     636                  spole = spole + varyear(i,nbp_lat,k,imth)
     637               END DO
     638               npole = npole/REAL(nbp_lon)
     639               spole = spole/REAL(nbp_lon)
     640               varyear(:,1,    k,imth) = npole
     641               varyear(:,nbp_lat,k,imth) = spole
     642            END DO
     643         END DO ! imth
    534644       
    535        ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
    536        IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
     645         ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
     646         IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
    537647       
    538        ! Transform from 2D to 1D field
    539        CALL grid2Dto1D_glo(varyear,varyear_glo1D)
    540        CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
    541        CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
    542 
     648         ! Transform from 2D to 1D field
     649         CALL grid2Dto1D_glo(varyear,varyear_glo1D)
     650         CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
     651         CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
     652     
     653      ENDIF
     654     
    543655    ELSE
    544       ALLOCATE(varyear_glo1D(0,0,0))       
     656        ALLOCATE(varyear_glo1D(0,0,0))       
    545657    END IF ! is_mpi_root .AND. is_omp_root
    546658
    … …  
    566678    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
    567679    ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
     680    ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr)
    568681    IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
    569682
    570     ! Scatter global field to local domain at local process
    571     CALL scatter(varyear_glo1D, pt_year)
    572     CALL scatter(psurf_glo1D, psurf_out)
    573     CALL scatter(load_glo1D,  load_out)
    574 
     683    IF (grid_type==unstructured) THEN
     684      IF (is_omp_master) THEN
     685        CALL xios_send_field(TRIM(varname)//"_in",varyear)
     686        CALL xios_recv_field(TRIM(varname)//"_out",pt_year_mpi)
     687        CALL xios_send_field("load_"//TRIM(varname)//"_in",load_glo2D)
     688        CALL xios_recv_field("load_"//TRIM(varname)//"_out",load_out_mpi)
     689        IF (first) THEN
     690          ALLOCATE(psurf_interp(klon_mpi,12))
     691          CALL xios_send_field("psurf_aerosol_in",psurf_glo2D)
     692          CALL xios_recv_field("psurf_aerosol_out",psurf_interp)
     693        ENDIF
     694      ENDIF
     695      CALL scatter_omp(pt_year_mpi,pt_year)
     696      CALL scatter_omp(load_out_mpi,load_out)
     697      CALL scatter_omp(psurf_interp,psurf_out)
     698      first=.FALSE.
     699    ELSE
     700      ! Scatter global field to local domain at local process
     701      CALL scatter(varyear_glo1D, pt_year)
     702      CALL scatter(psurf_glo1D, psurf_out)
     703      CALL scatter(load_glo1D,  load_out)
     704    ENDIF
    575705! 7) Test for negative values
    576706!****************************************************************************************

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosolstrato.F90

    r2745 r3413  
    77    USE phys_cal_mod, ONLY : mth_cur
    88    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, &
    9                                  grid2dto1d_glo
     9                                 grid2dto1d_glo, grid_type, unstructured
    1010    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
    1111    USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
    … …  
    1515    USE aero_mod
    1616    USE dimphy
     17    USE xios
    1718
    1819    implicit none
    … …  
    4344    real, allocatable:: tauaerstrat_mois(:, :, :)
    4445    real, allocatable:: tauaerstrat_mois_glo(:, :)
     46    real, allocatable:: tau_aer_strat_mpi(:, :)
    4547
    4648! For NetCDF:
    … …  
    8789    n_lat = size(latitude)
    8890    print *, 'LAT aerosol strato=', n_lat, latitude
    89     IF (n_lat.NE.nbp_lat) THEN
    90        print *,'Le nombre de lat n est pas egal a nbp_lat'
    91        STOP
    92     ENDIF
    93 
     91    IF (grid_type/=unstructured) THEN
     92      IF (n_lat.NE.nbp_lat) THEN
     93         print *,'Le nombre de lat n est pas egal a nbp_lat'
     94         STOP
     95      ENDIF
     96    ENDIF
     97   
    9498    CALL nf95_inq_varid(ncid_in, "LON", varid)
    9599    CALL nf95_gw_var(ncid_in, varid, longitude)
    96100    n_lon = size(longitude)
    97     print *, 'LON aerosol strato=', n_lon, longitude
    98     IF (n_lon.NE.nbp_lon) THEN
    99        print *,'Le nombre de lon n est pas egal a nbp_lon'
    100        STOP
    101     ENDIF
    102 
     101    IF (grid_type/=unstructured) THEN
     102      print *, 'LON aerosol strato=', n_lon, longitude
     103      IF (n_lon.NE.nbp_lon) THEN
     104         print *,'Le nombre de lon n est pas egal a nbp_lon'
     105         STOP
     106      ENDIF
     107    ENDIF
     108   
    103109    CALL nf95_inq_varid(ncid_in, "TIME", varid)
    104110    CALL nf95_gw_var(ncid_in, varid, time)
    … …  
    130136    CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
    131137
     138    ELSE
     139      ALLOCATE(tauaerstrat_mois(0,0,0))
    132140    ENDIF !--is_mpi_root and is_omp_root
    133141
    134142!$OMP BARRIER
    135143
     144    IF (grid_type==unstructured) THEN
     145      IF (is_omp_master) THEN
     146        CALL xios_send_field("taustrat_in",tauaerstrat_mois)
     147        ALLOCATE(tau_aer_strat_mpi(klon_mpi, klev))
     148        CALL xios_recv_field("taustrat_out",tau_aer_strat_mpi)
     149      ELSE
     150        ALLOCATE(tau_aer_strat_mpi(0,0))
     151      ENDIF
     152      CALL scatter_omp(tau_aer_strat_mpi,tau_aer_strat)
     153    ELSE 
    136154!--scatter on all proc
    137     CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
     155      CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
     156    ENDIF
    138157
    139158!--keep memory of previous month

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_horiz_time_climoz_m.F90

    r3411 r3413  
    22
    33  USE interpolation,     ONLY: locate
    4   USE mod_grid_phy_lmdz, ONLY: nlon_ou => nbp_lon, nlat_ou => nbp_lat
     4  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, grid_type, unstructured
    55  USE nrtype,            ONLY: pi
    66  USE netcdf,   ONLY: NF90_CLOBBER, NF90_FLOAT,     NF90_GET_VAR, NF90_OPEN,   &
    … …  
    1111          NF95_CLOSE, NF95_ENDDEF,  NF95_PUT_ATT,   NF95_PUT_VAR, NF95_COPY_ATT
    1212  USE print_control_mod, ONLY: lunout
     13  USE dimphy
    1314  IMPLICIT NONE
    1415  PRIVATE
    … …  
    1617  REAL, PARAMETER :: deg2rad=pi/180.
    1718  CHARACTER(LEN=13), PARAMETER :: vars_in(2)=['tro3         ','tro3_daylight']
     19
     20  INTEGER :: nlat_ou, nlon_ou
     21  REAL, ALLOCATABLE :: latitude_glo(:)
     22!$OMP THREADPRIVATE(latitude_glo)
     23  INTEGER, ALLOCATABLE :: ind_cell_glo_glo(:)
     24!$OMP THREADPRIVATE(ind_cell_glo_glo)
    1825
    1926CONTAINS
    … …  
    5259  USE assert_m,           ONLY: assert
    5360  USE cal_tools_m,        ONLY: year_len, mid_month
     61!!  USE control_mod,        ONLY: anneeref
    5462  USE time_phylmdz_mod,   ONLY: annee_ref
    5563  USE ioipsl,             ONLY: ioget_year_len, ioget_calendar
    … …  
    5866  USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
    5967  USE slopes_m,           ONLY: slopes
     68  USE xios
     69  USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
     70  USE geometry_mod, ONLY : latitude_deg, ind_cell_glo
     71  USE mod_grid_phy_lmdz, ONLY: klon_glo
     72
    6073!-------------------------------------------------------------------------------
    6174! Arguments:
    … …  
    8396  CHARACTER(LEN=20) :: cal_in              ! Calendar
    8497  REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:)   ! Ozone climatologies
     98  REAL, ALLOCATABLE :: o3_in3bis(:,:,:,:,:)   ! Ozone climatologies
    8599  REAL, ALLOCATABLE :: o3_in2  (:,:,:,:)   ! Ozone climatologies
     100  REAL, ALLOCATABLE :: o3_in2bis(:,:,:,:,:)   ! Ozone climatologies
    86101  ! last index: 1 for the day-night average, 2 for the daylight field.
    87102  REAL :: NaN
    … …  
    91106  REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13   ,:)
    92107  REAL, ALLOCATABLE :: o3_out3       (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
     108  REAL, ALLOCATABLE :: o3_out3_glo   (:,:,:,:) !   (nbp_lat,:,ntim_ou,:)
    93109  REAL, ALLOCATABLE :: o3_regr_lat     (:,:,:,:) !         (nlat_in,:,0:13   ,:)
    94110  REAL, ALLOCATABLE :: o3_out2         (:,:,:,:) !         (nlat_ou,:,ntim_ou,:)
     111  REAL, ALLOCATABLE :: o3_out2_glo     (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
     112  REAL, ALLOCATABLE :: o3_out          (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
    95113! Dimension number  | Interval                | Contains  | For variables:
    96114!   1 (longitude)   | [rlonu(i-1), rlonu(i)]  | rlonv(i)  | all
    … …  
    116134  INTEGER, ALLOCATABLE :: sta(:), cnt(:)
    117135  CHARACTER(LEN=80) :: sub, dim_nam, msg
    118 !-------------------------------------------------------------------------------
    119   sub="regr_horiz_time_climoz"
    120   WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
    121   CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
    122 
    123   CALL  NF95_OPEN("climoz.nc"  , NF90_NOWRITE, fID_in)
    124   lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
    125   lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
    126 
    127   !--- Get coordinates from the input file. Converts lon/lat in radians.
    128   !    Few inversions because "regr_conserv" and gcm need ascending vectors.
    129   CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
    130   CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
    131   l3D=ndims==4; l2D=ndims==3
    132   IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
    133   IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
    134   DO i=1,ndims
    135     CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
    136     CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
    137     ii=i; IF(l2D) ii=i+1                              !--- ndims==3:NO LONGITUDE
    138     SELECT CASE(ii)
    139       CASE(1)                                         !--- LONGITUDE
    140         CALL NF95_GW_VAR(fID_in, varid, lon_in)
    141         ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
    142         nlon_in=dln; lon_in=lon_in*deg2rad
    143       CASE(2)                                         !--- LATITUDE
    144         CALL NF95_GW_VAR(fID_in, varid, lat_in)
    145         ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
    146         nlat_in=dln; lat_in=lat_in*deg2rad
    147       CASE(3)                                         !--- PRESSURE LEVELS
    148         CALL NF95_GW_VAR(fID_in, varid, lev_in)
    149         ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
    150         nlev_in=dln
    151         CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
    152         k=LEN_TRIM(press_unit)
    153         DO WHILE(ICHAR(press_unit(k:k))==0)
    154           press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
    155         END DO
    156         IF(press_unit ==  "Pa") THEN
    157           lev_in = lev_in/100.                        !--- CONVERT TO hPa
    158         ELSE IF(press_unit /= "hPa") THEN
    159           CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
     136  REAL :: null_array(0)
     137  LOGICAL,SAVE :: first=.TRUE.
     138!$OMP THREADPRIVATE(first) 
     139  REAL, ALLOCATABLE :: test_o3_in(:,:)
     140  REAL, ALLOCATABLE :: test_o3_out(:)
     141
     142
     143  IF (grid_type==unstructured) THEN
     144    IF (first) THEN
     145      IF (is_master) THEN
     146        ALLOCATE(latitude_glo(klon_glo))
     147        ALLOCATE(ind_cell_glo_glo(klon_glo))
     148      ELSE
     149        ALLOCATE(latitude_glo(0))
     150        ALLOCATE(ind_cell_glo_glo(0))
     151      ENDIF
     152      CALL gather(latitude_deg,  latitude_glo)
     153      CALL gather(ind_cell_glo,  ind_cell_glo_glo)
     154    ENDIF
     155  ENDIF
     156   
     157  IF (is_omp_master) THEN
     158    nlat_ou=nbp_lat
     159    nlon_ou=nbp_lon
     160   
     161   !-------------------------------------------------------------------------------
     162    IF (is_mpi_root) THEN
     163      sub="regr_horiz_time_climoz"
     164      WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
     165      CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
     166
     167      CALL  NF95_OPEN("climoz.nc"  , NF90_NOWRITE, fID_in)
     168      lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
     169      lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
     170
     171      !--- Get coordinates from the input file. Converts lon/lat in radians.
     172      !    Few inversions because "regr_conserv" and gcm need ascending vectors.
     173      CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
     174      CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
     175      l3D=ndims==4; l2D=ndims==3
     176      IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
     177      IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
     178      DO i=1,ndims
     179        CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
     180        CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
     181        ii=i; IF(l2D) ii=i+1                              !--- ndims==3:NO LONGITUDE
     182        SELECT CASE(ii)
     183          CASE(1)                                         !--- LONGITUDE
     184            CALL NF95_GW_VAR(fID_in, varid, lon_in)
     185            ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
     186            nlon_in=dln; lon_in=lon_in*deg2rad
     187          CASE(2)                                         !--- LATITUDE
     188            CALL NF95_GW_VAR(fID_in, varid, lat_in)
     189            ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
     190            nlat_in=dln; lat_in=lat_in*deg2rad
     191          CASE(3)                                         !--- PRESSURE LEVELS
     192            CALL NF95_GW_VAR(fID_in, varid, lev_in)
     193            ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
     194            nlev_in=dln
     195            CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
     196            k=LEN_TRIM(press_unit)
     197            DO WHILE(ICHAR(press_unit(k:k))==0)
     198              press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
     199            END DO
     200            IF(press_unit ==  "Pa") THEN
     201              lev_in = lev_in/100.                        !--- CONVERT TO hPa
     202            ELSE IF(press_unit /= "hPa") THEN
     203              CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
     204            END IF
     205          CASE(4)                                         !--- TIME
     206            CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
     207            cal_in='gregorian'
     208            IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR)        & 
     209            WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'//       &
     210            TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
     211            k=LEN_TRIM(cal_in)
     212            DO WHILE(ICHAR(cal_in(k:k))==0)
     213              cal_in(k:k)=' '; k=LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
     214            END DO
     215        END SELECT
     216      END DO
     217
     218      !--- Prepare quantities for time interpolation
     219      tmidmonth=mid_month(annee_ref, cal_in)
     220      IF(interpt) THEN
     221        ntim_ou=ioget_year_len(annee_ref)
     222        ALLOCATE(tmidday(ntim_ou))
     223        tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
     224        CALL ioget_calendar(cal_ou)
     225      ELSE
     226        ntim_ou=14
     227        cal_ou=cal_in
     228      END IF
     229    ENDIF
     230
     231    IF (grid_type==unstructured) THEN
     232      CALL bcast_mpi(nlon_in)
     233      CALL bcast_mpi(nlat_in)
     234      CALL bcast_mpi(nlev_in)
     235      CALL bcast_mpi(l3d)
     236      CALL bcast_mpi(tmidmonth)
     237      CALL bcast_mpi(tmidday)
     238      CALL bcast_mpi(ntim_ou)
     239   
     240      IF (is_mpi_root) THEN
     241        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=nlat_in, jbegin=0, latvalue_1d=lat_in/deg2rad)
     242        IF (l3D) THEN
     243          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=nlon_in, ibegin=0, lonvalue_1d=lon_in/deg2rad)
     244        ELSE
     245          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=8, ibegin=0, lonvalue_1d = (/ 0.,45.,90.,135.,180.,225.,270., 315. /))
     246        ENDIF
     247      ELSE
     248        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=0, jbegin=0, latvalue_1d=null_array )
     249        IF (l3D) THEN
     250          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=0, ibegin=0, lonvalue_1d=null_array)
     251        ELSE
     252          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=0, ibegin=0, lonvalue_1d=null_array)
     253        ENDIF
     254      ENDIF
     255      CALL  xios_set_axis_attr("axis_climoz", n_glo=nlev_in)
     256      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
     257      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
     258      CALL  xios_set_axis_attr("tr_climoz", n_glo=read_climoz)
     259      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
     260      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
     261
     262      IF (first) THEN
     263        first=.FALSE.
     264        RETURN
     265      ENDIF
     266    ENDIF
     267   
     268   
     269    IF (is_mpi_root) THEN     
     270      !--- Longitudes management:
     271      !    * Need to shift data if the origin of input file longitudes /= -pi
     272      !    * Need to add some margin in longitude to ensure input interval contains
     273      !      all the output intervals => at least one longitudes slice has to be
     274      !      duplicated, possibly more for undersampling.
     275      IF(l3D) THEN
     276        IF (grid_type==unstructured) THEN
     277          dx2=0
     278        ELSE
     279          !--- Compute input edges longitudes vector (no end point yet)
     280          ALLOCATE(v1(nlon_in+1))
     281          v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
     282          FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
     283          v1(nlon_in+1)=v1(1)+2.*pi
     284          DEALLOCATE(lon_in)
     285
     286          !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
     287          v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
     288
     289          !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
     290          dx1=locate(v1,boundslon_reg(1,west))-1
     291          v1=CSHIFT(v1,SHIFT=dx1,DIM=1)
     292          v1(nlon_in-dx1+1:)=v1(nlon_in-dx1+1:)+2.*pi
     293   
     294          !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
     295          dx2=0; DO WHILE(v1(1+dx2)+2.*pi<boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
     296
     297          !--- Final edges longitudes vector (with margin and end point)
     298          ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
     299          DEALLOCATE(v1)
     300        ENDIF
     301      END IF
     302
     303      !--- Compute sinus of intervals edges latitudes:
     304      ALLOCATE(sinlat_in_edge(nlat_in+1))
     305      sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
     306      FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
     307      DEALLOCATE(lat_in)
     308
     309
     310
     311      !--- Check for contiguous years:
     312      ib=0; ie=13
     313      IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
     314        WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
     315      ELSE 
     316        IF(     lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
     317        IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
     318        IF(     lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
     319        IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
     320        IF(.NOT.lprev) ib=1
     321        IF(.NOT.lnext) ie=12
     322      END IF
     323      ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1 
     324      IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
     325      IF(l2D) cnt=[        nlat_in,nlev_in,1] 
     326      IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
     327      IF(l2D) ALLOCATE(o3_in2(            nlat_in,nlev_in,ib:ie,read_climoz))
     328
     329      !--- Read full current file and one record each available contiguous file
     330      DO iv=1,read_climoz
     331        msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
     332        CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
     333        IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
     334        IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
     335        CALL handle_err(TRIM(msg), ncerr, fID_in)
     336        IF(lprev) THEN; sta(ndims)=12 
     337          CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
     338          IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
     339          IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
     340          CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
    160341        END IF
    161       CASE(4)                                         !--- TIME
    162         CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
    163         cal_in='gregorian'
    164         IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR)        &
    165           WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'//       &
    166           TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
    167         k=LEN_TRIM(cal_in)
    168         DO WHILE(ICHAR(cal_in(k:k))==0)
    169           cal_in(k:k)=' '; k=LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
    170         END DO
    171     END SELECT
    172   END DO
    173 
    174   !--- Longitudes management:
    175   !    * Need to shift data if the origin of input file longitudes /= -pi
    176   !    * Need to add some margin in longitude to ensure input interval contains
    177   !      all the output intervals => at least one longitudes slice has to be
    178   !      duplicated, possibly more for undersampling.
    179   IF(l3D) THEN
    180     !--- Compute input edges longitudes vector (no end point yet)
    181     ALLOCATE(v1(nlon_in+1))
    182     v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
    183     FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
    184     v1(nlon_in+1)=v1(1)+2.*pi
    185     DEALLOCATE(lon_in)
    186 
    187     !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
    188     v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
    189 
    190     !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
    191     dx1=locate(v1,boundslon_reg(1,west))-1
    192     v1=CSHIFT(v1,SHIFT=dx1,DIM=1); v1(nlon_in-dx1+1:)=v1(nlon_in-dx1+1:)+2.*pi
    193 
    194     !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
    195     dx2=0; DO WHILE(v1(1+dx2)+2.*pi<boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
    196 
    197     !--- Final edges longitudes vector (with margin and end point)
    198     ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
    199     DEALLOCATE(v1)
    200   END IF
    201 
    202   !--- Compute sinus of intervals edges latitudes:
    203   ALLOCATE(sinlat_in_edge(nlat_in+1))
    204   sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
    205   FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
    206   DEALLOCATE(lat_in)
    207 
    208   !--- Prepare quantities for time interpolation
    209   tmidmonth=mid_month(annee_ref, cal_in)
    210   IF(interpt) THEN
    211     ntim_ou=ioget_year_len(annee_ref)
    212     ALLOCATE(tmidday(ntim_ou))
    213     tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
    214     CALL ioget_calendar(cal_ou)
    215   ELSE
    216     ntim_ou=14
    217     cal_ou=cal_in
    218   END IF
    219 
    220   !--- Create the output file and get the variable IDs:
    221   CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
    222                    ndims, cal_ou)
    223 
    224   !--- Write remaining coordinate variables:
    225   CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
    226   IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
    227   IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
    228 
    229   !--- Check for contiguous years:
    230   ib=0; ie=13
    231   IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
    232     WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
    233   ELSE
    234     IF(     lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
    235     IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
    236     IF(     lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
    237     IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
    238     IF(.NOT.lprev) ib=1
    239     IF(.NOT.lnext) ie=12
    240   END IF
    241   ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1
    242   IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
    243   IF(l2D) cnt=[        nlat_in,nlev_in,1]
    244   IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
    245   IF(l2D) ALLOCATE(o3_in2(            nlat_in,nlev_in,ib:ie,read_climoz))
    246 
    247   !--- Read full current file and one record each available contiguous file
    248   DO iv=1,read_climoz
    249     msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
    250     CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
    251     IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
    252     IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
    253     CALL handle_err(TRIM(msg), ncerr, fID_in)
    254     IF(lprev) THEN; sta(ndims)=12
    255       CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
    256       IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
    257       IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
    258       CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
     342        IF(lnext) THEN; sta(ndims)=1 
     343          CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
     344          IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
     345          IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
     346          CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
     347        END IF
     348      END DO
     349      IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
     350      IF(lprev) CALL NF95_CLOSE(fID_in_m)
     351      IF(lnext) CALL NF95_CLOSE(fID_in_p)
     352
     353      !--- Revert decreasing coordinates vector
     354      IF(l3D) THEN
     355        IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
     356        IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
     357        IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
     358       
     359        IF (grid_type /= unstructured) THEN
     360          !--- Shift values for longitude and duplicate some longitudes slices
     361          o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
     362          o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
     363        ENDIF
     364      ELSE
     365        IF(ldec_lat) o3_in2 = o3_in2(  nlat_in:1:-1,:,:,:)
     366        IF(ldec_lev) o3_in2 = o3_in2(  :,nlev_in:1:-1,:,:)
     367      END IF
     368
     369     !--- Deal with missing values
     370      DO m=1, read_climoz
     371        WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
     372        IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
     373          IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
     374            WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
     375          END IF
     376        END IF
     377        WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
     378        WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
     379
     380        !--- Check top layer contains no NaNs & search NaNs from top to ground
     381        msg=TRIM(sub)//": NaNs in top layer !"
     382        IF(l3D) THEN
     383          IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
     384          DO k = 2,nlev_in
     385            WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
     386          END DO
     387        ELSE
     388          IF(ANY(o3_in2(  :,1,:,m)==NaN)) THEN
     389            WRITE(lunout,*)msg
     390            !--- Fill in latitudes where all values are missing
     391            DO l=1,nmth_in
     392              !--- Next to south pole
     393              j=1;       DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
     394              IF(j>1) &
     395                o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
     396              !--- Next to north pole
     397              j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
     398              IF(j<nlat_in) &
     399                o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
     400            END DO
     401          END IF
     402
     403          !--- Fill in high latitudes missing values
     404          !--- Highest level been filled-in, so has always valid values.
     405          DO k = 2,nlev_in
     406            WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
     407          END DO
     408        END IF
     409      END DO
     410
     411    ENDIF
     412   
     413    !=============================================================================
     414    IF(l3D) THEN                                                   !=== 3D FIELDS
     415    !=============================================================================
     416     IF (grid_type==unstructured) THEN
     417       nlat_ou=klon_mpi
     418       
     419       IF (is_mpi_root) THEN
     420         ALLOCATE(o3_in3bis(nlon_in,nlat_in,nlev_in,0:13,read_climoz))
     421         o3_in3bis(:,:,:,ib:ie,:)=o3_in3(1:nlon_in,:,:,ib:ie,:)
     422       ELSE
     423         ALLOCATE(o3_in3bis(0,0,0,0,read_climoz))
     424       ENDIF
     425       ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
     426       
     427       CALL xios_send_field("tro3_in",o3_in3bis(:,:,:,:,:))
     428       CALL xios_recv_field("tro3_out",o3_regr_lonlat(1,:,:,:,:))
     429
     430     ELSE
     431         
     432       !--- Regrid in longitude
     433        ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
     434        CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
     435                            xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
     436                            vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
     437        DEALLOCATE(o3_in3)
     438
     439        !--- Regrid in latitude: averaging with respect to SIN(lat) is
     440        !                        equivalent to weighting by COS(lat)
     441        !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
     442        ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
     443        CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
     444                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
     445                        vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
     446                   slope = slopes(2,o3_regr_lon, sinlat_in_edge))
     447        DEALLOCATE(o3_regr_lon)
     448
     449     ENDIF
     450
     451     !--- Duplicate previous/next record(s) if they are not available
     452     IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
     453     IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
     454     
     455     !--- Regrid in time by linear interpolation:
     456     ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
     457     IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
     458     IF(.NOT.interpt) o3_out3=o3_regr_lonlat
     459     DEALLOCATE(o3_regr_lonlat)
     460
     461     nlat_ou=nbp_lat
     462     IF (grid_type==unstructured) THEN
     463       CALL xios_send_field('o3_out',o3_out3)
     464       ndims=3
     465       ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
     466       CALL gather_mpi(o3_out3(1,:,:,:,:), o3_out3_glo)
     467     ENDIF
     468
     469    !--- Create the output file and get the variable IDs:
     470    CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
     471                     ndims, cal_ou)
     472
     473    IF (is_mpi_root) THEN
     474      !--- Write remaining coordinate variables:
     475      CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
     476      IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
     477      IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
     478
     479      !--- Write to file (the order of "rlatu" is inverted in the output file):
     480        IF (grid_type==unstructured) THEN
     481
     482          ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
     483          DO i=1,klon_glo
     484            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out3_glo(i,:,:,:)
     485          ENDDO
     486
     487          DO m = 1, read_climoz
     488            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
     489          END DO
     490         
     491        ELSE
     492          DO m = 1, read_climoz
     493            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
     494          END DO
     495      ENDIF
     496      CALL NF95_CLOSE(fID_ou)
     497
     498
     499    ENDIF
     500
     501
     502    !=============================================================================
     503    ELSE                                                         !=== ZONAL FIELDS
     504    !=============================================================================
     505   
     506     IF (grid_type==unstructured) THEN
     507       nlat_ou=klon_mpi
     508
     509       IF (is_mpi_root) THEN
     510         ALLOCATE(o3_in2bis(8,nlat_in,nlev_in,0:13,read_climoz))
     511         o3_in2bis(:,:,:,ib:ie,:)=SPREAD(o3_in2,1,8)
     512       ELSE
     513         ALLOCATE(o3_in2bis(0,0,0,0,read_climoz))
     514       ENDIF
     515       ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
     516       CALL xios_send_field("tro3_in",o3_in2bis(:,:,:,:,:))
     517       CALL xios_recv_field("tro3_out",o3_regr_lat(:,:,:,:))
     518       IF(.NOT.lprev) o3_regr_lat(:,:, 0, :) = o3_regr_lat(:,:,12,:)
     519       IF(.NOT.lnext) o3_regr_lat(:,:,13, :) = o3_regr_lat(:,:, 1,:)
     520       
     521     
     522     ELSE
     523        !--- Regrid in latitude: averaging with respect to SIN(lat) is
     524        !                        equivalent to weighting by COS(lat)
     525        !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
     526        ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
     527        CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
     528                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
     529                        vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
     530                     slope = slopes(1,o3_in2, sinlat_in_edge))
     531        DEALLOCATE(o3_in2)
     532
     533        !--- Duplicate previous/next record(s) if they are not available
     534        IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
     535        IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
     536
     537     ENDIF
     538     
     539      !--- Regrid in time by linear interpolation:
     540      ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
     541      IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
     542      IF(.NOT.interpt) o3_out2=o3_regr_lat
     543      DEALLOCATE(o3_regr_lat)
     544
     545      nlat_ou=nbp_lat
     546   
     547      IF (grid_type==unstructured) THEN
     548        ndims=3
     549        ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
     550        CALL gather_mpi(o3_out2, o3_out2_glo)
     551      ENDIF
     552     
     553      !--- Create the output file and get the variable IDs:
     554      CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
     555                         ndims, cal_ou)
     556
     557      IF (is_mpi_root) THEN
     558     
     559        !--- Write remaining coordinate variables:
     560        CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
     561        IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
     562        IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
     563
     564        IF (grid_type==unstructured) THEN
     565
     566          ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
     567          DO i=1,klon_glo
     568            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out2_glo(i,:,:,:)
     569          ENDDO
     570
     571
     572          DO m = 1, read_climoz
     573            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
     574          END DO
     575        ELSE
     576          !--- Write to file (the order of "rlatu" is inverted in the output file):
     577          DO m = 1, read_climoz
     578            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
     579          END DO
     580        ENDIF
     581       
     582        CALL NF95_CLOSE(fID_ou)
     583     
     584      ENDIF
     585
     586    !=============================================================================
    259587    END IF
    260     IF(lnext) THEN; sta(ndims)=1
    261       CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
    262       IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
    263       IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
    264       CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
    265     END IF
    266   END DO
    267   IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
    268   IF(lprev) CALL NF95_CLOSE(fID_in_m)
    269   IF(lnext) CALL NF95_CLOSE(fID_in_p)
    270 
    271   !--- Revert decreasing coordinates vector
    272   IF(l3D) THEN
    273     IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
    274     IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
    275     IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
    276     !--- Shift values for longitude and duplicate some longitudes slices
    277     o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
    278     o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
    279   ELSE
    280     IF(ldec_lat) o3_in2 = o3_in2(  nlat_in:1:-1,:,:,:)
    281     IF(ldec_lev) o3_in2 = o3_in2(  :,nlev_in:1:-1,:,:)
    282   END IF
    283 
    284  !--- Deal with missing values
    285   DO m=1, read_climoz
    286     WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
    287     IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
    288       IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
    289         WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
    290       END IF
    291     END IF
    292     WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
    293     WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
    294 
    295     !--- Check top layer contains no NaNs & search NaNs from top to ground
    296     msg=TRIM(sub)//": NaNs in top layer !"
    297     IF(l3D) THEN
    298       IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
    299       DO k = 2,nlev_in
    300         WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
    301       END DO
    302     ELSE
    303       IF(ANY(o3_in2(  :,1,:,m)==NaN)) THEN
    304         WRITE(lunout,*)msg
    305         !--- Fill in latitudes where all values are missing
    306         DO l=1,nmth_in
    307           !--- Next to south pole
    308           j=1;       DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
    309           IF(j>1) &
    310             o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
    311           !--- Next to north pole
    312           j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
    313           IF(j<nlat_in) &
    314             o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
    315         END DO
    316       END IF
    317 
    318       !--- Fill in high latitudes missing values
    319       !--- Highest level been filled-in, so has always valid values.
    320       DO k = 2,nlev_in
    321         WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
    322       END DO
    323     END IF
    324   END DO
    325   CALL NF95_CLOSE(fID_in)
    326 
    327   !=============================================================================
    328   IF(l3D) THEN                                                   !=== 3D FIELDS
    329   !=============================================================================
    330     !--- Regrid in longitude
    331     ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
    332     CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
    333                         xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
    334                         vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
    335     DEALLOCATE(o3_in3)
    336 
    337     !--- Regrid in latitude: averaging with respect to SIN(lat) is
    338     !                        equivalent to weighting by COS(lat)
    339     !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
    340     ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
    341     CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
    342                     xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
    343                     vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
    344                  slope = slopes(2,o3_regr_lon, sinlat_in_edge))
    345     DEALLOCATE(o3_regr_lon)
    346 
    347     !--- Duplicate previous/next record(s) if they are not available
    348     IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
    349     IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
    350 
    351     !--- Regrid in time by linear interpolation:
    352     ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
    353     IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
    354     IF(.NOT.interpt) o3_out3=o3_regr_lonlat
    355     DEALLOCATE(o3_regr_lonlat)
    356 
    357     !--- Write to file (the order of "rlatu" is inverted in the output file):
    358     DO m = 1, read_climoz
    359       CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
    360     END DO
    361 
    362   !=============================================================================
    363   ELSE                                                         !=== ZONAL FIELDS
    364   !=============================================================================
    365     !--- Regrid in latitude: averaging with respect to SIN(lat) is
    366     !                        equivalent to weighting by COS(lat)
    367     !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
    368     ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
    369     CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
    370                     xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
    371                     vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
    372                  slope = slopes(1,o3_in2, sinlat_in_edge))
    373     DEALLOCATE(o3_in2)
    374 
    375     !--- Duplicate previous/next record(s) if they are not available
    376     IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
    377     IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
    378 
    379     !--- Regrid in time by linear interpolation:
    380     ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
    381     IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
    382     IF(.NOT.interpt) o3_out2=o3_regr_lat
    383     DEALLOCATE(o3_regr_lat)
    384 
    385     !--- Write to file (the order of "rlatu" is inverted in the output file):
    386     DO m = 1, read_climoz
    387       CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
    388     END DO
    389 
    390   !=============================================================================
    391   END IF
    392   !=============================================================================
    393 
    394   CALL NF95_CLOSE(fID_ou)
    395 
     588    !=============================================================================
     589
     590    IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
     591
     592  ENDIF ! is_omp_master
     593
     594  first=.FALSE.
    396595END SUBROUTINE regr_horiz_time_climoz
    397596!
    … …  
    408607!-------------------------------------------------------------------------------
    409608  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
     609  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
     610  USE mod_phys_lmdz_para, ONLY: is_mpi_root
     611  USE mod_grid_phy_lmdz, ONLY: klon_glo
     612!
    410613!-------------------------------------------------------------------------------
    411614! Arguments:
    … …  
    419622  INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
    420623  INTEGER :: vlonID, vlatID, ncerr,  is
     624  REAL,ALLOCATABLE    :: latitude_glo_(:)
    421625  CHARACTER(LEN=80) :: sub
    422 !-------------------------------------------------------------------------------
    423   sub="prepare_out"
    424   WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
    425   CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
     626  INTEGER :: i
     627
     628
     629!-------------------------------------------------------------------------------
     630 
     631  IF (is_mpi_root) THEN 
     632    sub="prepare_out"
     633    WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
     634    CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
    426635
    427636  !--- Dimensions:
    428   IF(ndims==4) &
    429   CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
    430   CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
    431   CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
    432   CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
    433 
    434   !--- Define coordinate variables:
    435   IF(ndims==4) &
    436   CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
    437   CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
    438   CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
    439   CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
    440   IF(ndims==4) &
    441   CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
    442   CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
    443   CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
    444   CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
    445   IF(ndims==4) &
    446   CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
    447   CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
    448   CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
    449   CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
    450   CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
    451   CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
     637    IF(ndims==4) &
     638    CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
     639    CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
     640    CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
     641    CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
     642
     643    !--- Define coordinate variables:
     644    IF(ndims==4) &
     645    CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
     646    CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
     647    CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
     648    CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
     649    IF(ndims==4) &
     650    CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
     651    CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
     652    CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
     653    CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
     654    IF(ndims==4) &
     655    CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
     656    CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
     657    CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
     658    CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
     659    CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
     660    CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
    452661
    453662  !--- Define the main variables:
    454   IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
    455   IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
    456   CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
    457   CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
    458   CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
     663    IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
     664    IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
     665    CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
     666    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
     667    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
    459668      &_in_air")
    460   IF(SIZE(vID_ou) == 2) THEN
    461     CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
    462     CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
    463       &ylight")
    464   END IF
     669    IF(SIZE(vID_ou) == 2) THEN
     670      CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
     671      CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
     672        &ylight")
     673    END IF
    465674
    466675  !--- Global attributes:
    467676  ! The following commands, copying attributes, may fail. That is OK.
    468677  ! It should just mean that the attribute is not defined in the input file.
    469   CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
    470   CALL handle_err_copy_att("Conventions")
    471   CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
    472   CALL handle_err_copy_att("title")
    473   CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
    474   CALL handle_err_copy_att("institution")
    475   CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
    476   CALL handle_err_copy_att("source")
    477   CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
    478   CALL NF95_ENDDEF(fID_ou)
    479 
    480   !--- Write one of the coordinate variables:
    481   IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
    482   CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
    483   !    (convert from rad to degrees and sort in ascending order)
    484 
     678    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
     679    CALL handle_err_copy_att("Conventions")
     680    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
     681    CALL handle_err_copy_att("title")
     682    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
     683    CALL handle_err_copy_att("institution")
     684    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
     685    CALL handle_err_copy_att("source")
     686    CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
     687    CALL NF95_ENDDEF(fID_ou)
     688
     689    IF (grid_type==unstructured) THEN
     690      ALLOCATE(latitude_glo_(klon_glo))
     691      DO i=1,klon_glo
     692        latitude_glo_(ind_cell_glo_glo(i))=latitude_glo(i)
     693      ENDDO
     694      CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
     695    ELSE
     696      !--- Write one of the coordinate variables:
     697      IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
     698      CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
     699    !    (convert from rad to degrees and sort in ascending order)
     700    ENDIF
     701  ENDIF
     702 
    485703CONTAINS
    486704

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_pr_time_av_m.F90

    r3411 r3413  
    106106  USE assert_m,       ONLY: assert
    107107  USE assert_eq_m,    ONLY: assert_eq
    108 !  USE comvert_mod,    ONLY: scaleheight
     108!!  USE comvert_mod,    ONLY: scaleheight
    109109  USE interpolation,  ONLY: locate
    110110  USE regr_conserv_m, ONLY: regr_conserv
    111111  USE regr_lint_m,    ONLY: regr_lint
    112112  USE slopes_m,       ONLY: slopes
    113   USE mod_phys_lmdz_mpi_data,       ONLY: is_mpi_root
    114   USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev
    115   USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter
     113  USE mod_phys_lmdz_para,           ONLY: is_mpi_root,is_master
     114  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid_type, unstructured
     115  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter, gather
    116116  USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
     117  USE geometry_mod,                 ONLY: ind_cell_glo
    117118!-------------------------------------------------------------------------------
    118119! Arguments:
    … …  
    157158  REAL, DIMENSION(nbp_lon, nbp_lat)   :: ps1, pt1, ot1
    158159  REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
     160  INTEGER,ALLOCATABLE :: ind_cell_glo_glo(:)
    159161  LOGICAL :: ll
    160162!-------------------------------------------------------------------------------
    … …  
    230232    CALL bcast(lO3Trop); CALL bcast(linterp)
    231233  END IF
     234 
     235  IF (is_master) THEN
     236    ALLOCATE(ind_cell_glo_glo(klon_glo))
     237  ELSE
     238    ALLOCATE(ind_cell_glo_glo(0))
     239  ENDIF
     240  CALL gather(ind_cell_glo,ind_cell_glo_glo)
     241  IF (is_master .AND. grid_type==unstructured) v1(:,:,:,:)=v1(:,ind_cell_glo_glo(:),:,:)
     242 
    232243  CALL scatter2d(v1,v2)
     244
    233245  !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
    234   IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF
    235   IF(lPrTrop) CALL scatter2d(pt1,pt2)
    236   IF(lO3Trop) CALL scatter2d(ot1,ot2)
    237 
     246  IF(lPrSurf) THEN
     247    IF (is_master .AND. grid_type==unstructured) ps1(:,:)=ps1(:,ind_cell_glo_glo(:))
     248    CALL scatter2d(ps1,ps2)
     249  ELSE
     250   ps2=pint_ou(:,1)
     251  END IF
     252
     253  IF(lPrTrop) THEN
     254    IF (is_master .AND. grid_type==unstructured) pt1(:,:)=pt1(:,ind_cell_glo_glo(:))
     255    CALL scatter2d(pt1,pt2)
     256  ENDIF
     257
     258  IF(lO3Trop) THEN
     259    IF (is_master .AND. grid_type==unstructured) ot1(:,:)=ot1(:,ind_cell_glo_glo(:))
     260    CALL scatter2d(ot1,ot2)
     261  ENDIF
    238262  !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
    239263  IF(.NOT.lAdjTro) THEN

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato1_rrtm.F90

    r2744 r3413  
    22! $Id: readaerosolstrato1_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
    33!
     4
    45SUBROUTINE readaerosolstrato1_rrtm(debut)
    56
    … …  
    910
    1011    USE phys_cal_mod, ONLY : mth_cur
    11     USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo
    12     USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
    13     USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
     12    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo, grid_type, unstructured
    1413    USE mod_phys_lmdz_para
    1514    USE phys_state_var_mod
    … …  
    1918    USE YOERAD, ONLY : NLW
    2019    USE YOMCST
     20    USE xios
    2121
    2222    IMPLICIT NONE
    … …  
    4545    REAL, ALLOCATABLE:: tauaerstrat_mois(:, :, :)
    4646    REAL, ALLOCATABLE:: tauaerstrat_mois_glo(:, :)
     47    REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :)
    4748
    4849! For NetCDF:
    … …  
    102103    n_lat = size(latitude)
    103104    print *, 'LAT aerosol strato=', n_lat, latitude
    104     IF (n_lat.NE.nbp_lat) THEN
    105        print *,'Le nombre de lat n est pas egal a nbp_lat'
    106        STOP
    107     ENDIF
    108 
     105
     106    IF (grid_type/=unstructured) THEN
     107      IF (n_lat.NE.nbp_lat) THEN
     108         print *,'Le nombre de lat n est pas egal a nbp_lat'
     109         STOP
     110      ENDIF
     111    ENDIF
     112   
    109113    CALL nf95_inq_varid(ncid_in, "LON", varid)
    110114    CALL nf95_gw_var(ncid_in, varid, longitude)
    111115    n_lon = size(longitude)
    112116    print *, 'LON aerosol strato=', n_lon, longitude
    113     IF (n_lon.NE.nbp_lon) THEN
    114        print *,'Le nombre de lon n est pas egal a nbp_lon'
    115        STOP
    116     ENDIF
    117 
     117
     118    IF (grid_type/=unstructured) THEN
     119      IF (n_lon.NE.nbp_lon) THEN
     120         print *,'Le nombre de lon n est pas egal a nbp_lon'
     121         STOP
     122      ENDIF
     123    ENDIF
     124   
     125   
    118126    CALL nf95_inq_varid(ncid_in, "TIME", varid)
    119127    CALL nf95_gw_var(ncid_in, varid, time)
    … …  
    144152!---reduce to a klon_glo grid
    145153    CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
    146 
     154   
     155    ELSE
     156      ALLOCATE(tauaerstrat_mois(0,0,0))
    147157    ENDIF !--is_mpi_root and is_omp_root
    148158
    … …  
    153163
    154164!--scatter on all proc
    155     CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
    156 
     165   
     166    IF (grid_type==unstructured) THEN
     167      IF (is_omp_master) THEN
     168        ALLOCATE(tauaerstrat_mpi(klon_mpi,klev))
     169        CALL xios_send_field("taustrat_in",tauaerstrat_mois)
     170        CALL xios_recv_field("taustrat_out",tauaerstrat_mpi)
     171      ELSE
     172        ALLOCATE(tauaerstrat_mpi(0,0))
     173      ENDIF
     174      CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
     175    ELSE 
     176      CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
     177    ENDIF
     178   
    157179    IF (is_mpi_root.AND.is_omp_root) THEN
    158180!

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90

    r2744 r3413  
    22! $Id: readaerosolstrato2_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
    33!
     4
    45SUBROUTINE readaerosolstrato2_rrtm(debut)
    56
    … …  
    910
    1011    USE phys_cal_mod, ONLY : mth_cur
    11     USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo
    12     USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
    13     USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
     12    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo, grid_type, unstructured
     13    USE mod_phys_lmdz_mpi_data
     14    USE mod_phys_lmdz_omp_data
    1415    USE mod_phys_lmdz_para
    1516    USE phys_state_var_mod
    … …  
    1920    USE YOERAD, ONLY : NLW
    2021    USE YOMCST
    21 
     22    USE xios
    2223    IMPLICIT NONE
    2324
    … …  
    6566    REAL, ALLOCATABLE:: cgaerstrat_mois_glo(:, :, :)
    6667    REAL, ALLOCATABLE:: taulwaerstrat_mois_glo(:, :, :)
     68    REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :, :)
     69    REAL, ALLOCATABLE:: pizaerstrat_mpi(:, :, :)
     70    REAL, ALLOCATABLE:: cgaerstrat_mpi(:, :, :)
     71    REAL, ALLOCATABLE:: taulwaerstrat_mpi(:, :, :)
    6772
    6873! For NetCDF:
    … …  
    107112        CALL nf95_gw_var(ncid_in, varid, latitude)
    108113        n_lat = size(latitude)
    109         IF (n_lat.NE.nbp_lat) THEN
    110            print *, 'latitude=', n_lat, nbp_lat
    111            abort_message='Le nombre de lat n est pas egal a nbp_lat'
    112            CALL abort_physic(modname,abort_message,1)
     114
     115        IF (grid_type/=unstructured) THEN
     116           IF (n_lat.NE.nbp_lat) THEN
     117             print *, 'latitude=', n_lat, nbp_lat
     118             abort_message='Le nombre de lat n est pas egal a nbp_lat'
     119             CALL abort_physic(modname,abort_message,1)
     120           ENDIF
    113121        ENDIF
    114122
    … …  
    134142        ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
    135143
    136         ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
    137         ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
    138         ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
    139 
    140         ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    141         ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    142         ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    143 
    144144!--reading stratospheric aerosol tau per layer
    145145        CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
    … …  
    159159        CALL nf95_close(ncid_in)
    160160
     161       
     162        IF (grid_type/=unstructured) THEN
     163          ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
     164          ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
     165          ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
     166
     167          ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
     168          ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
     169          ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    161170!--select the correct month
    162171!--and copy into 1st longitude
    163         tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
    164         pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
    165         cgaerstrat_mois(1,:,:,:)  = cgaerstrat(:,:,:,mth_cur)
     172          tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
     173          pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
     174          cgaerstrat_mois(1,:,:,:)  = cgaerstrat(:,:,:,mth_cur)
    166175
    167176!--copy longitudes
    168         DO i=2, n_lon
    169          tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
    170          pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
    171          cgaerstrat_mois(i,:,:,:)  = cgaerstrat_mois(1,:,:,:)
    172         ENDDO
     177          DO i=2, n_lon
     178           tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
     179           pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
     180           cgaerstrat_mois(i,:,:,:)  = cgaerstrat_mois(1,:,:,:)
     181          ENDDO
    173182
    174183!---reduce to a klon_glo grid
    175         DO band=1, NSW
    176           CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
    177           CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
    178           CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
    179         ENDDO
    180 
     184          DO band=1, NSW
     185            CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
     186            CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
     187            CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
     188          ENDDO
     189        ENDIF
    181190!--Now LW optical properties
    182191!
     192
    183193        CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
    184194
    … …  
    194204        CALL nf95_gw_var(ncid_in, varid, latitude)
    195205        n_lat = size(latitude)
    196         IF (n_lat.NE.nbp_lat) THEN
    197            abort_message='Le nombre de lat n est pas egal a nbp_lat'
    198            CALL abort_physic(modname,abort_message,1)
    199         ENDIF
    200 
     206
     207        IF (grid_type/=unstructured) THEN
     208          IF (n_lat.NE.nbp_lat) THEN
     209             abort_message='Le nombre de lat n est pas egal a nbp_lat'
     210             CALL abort_physic(modname,abort_message,1)
     211          ENDIF
     212        ENDIF
     213       
    201214        CALL nf95_inq_varid(ncid_in, "TIME", varid)
    202215        CALL nf95_gw_var(ncid_in, varid, time)
    … …  
    217230
    218231        ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
    219         ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
    220         ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
    221232
    222233!--reading stratospheric aerosol lw tau per layer
    … …  
    227238        CALL nf95_close(ncid_in)
    228239
     240        IF (grid_type/=unstructured) THEN
     241
     242          ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
     243          ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
     244
    229245!--select the correct month
    230246!--and copy into 1st longitude
    231         taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
     247          taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
    232248!--copy longitudes
    233         DO i=2, n_lon
    234           taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
    235         ENDDO
     249          DO i=2, n_lon
     250            taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
     251          ENDDO
    236252
    237253!---reduce to a klon_glo grid
    238         DO band=1, NLW
    239           CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
    240         ENDDO
    241 
     254          DO band=1, NLW
     255            CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
     256          ENDDO
     257        ENDIF
     258     
    242259      ELSE !--proc other than mpi_root and omp_root
    243260           !--dummy allocation needed for debug mode
    … …  
    248265        ALLOCATE(taulwaerstrat_mois_glo(1,1,1))
    249266
     267        ALLOCATE(tauaerstrat(0,0,0,12))
     268        ALLOCATE(pizaerstrat(0,0,0,12))
     269        ALLOCATE(cgaerstrat(0,0,0,12))
     270        ALLOCATE(taulwaerstrat(0,0,0,12))
     271
     272
    250273      ENDIF !--is_mpi_root and is_omp_root
    251274
    … …  
    255278      mth_pre=mth_cur
    256279
     280      IF (grid_type==unstructured) THEN
     281
     282        IF (is_omp_master) THEN
     283          ALLOCATE(tauaerstrat_mpi(klon_mpi, klev, NSW))
     284          ALLOCATE(pizaerstrat_mpi(klon_mpi, klev, NSW))
     285          ALLOCATE(cgaerstrat_mpi(klon_mpi, klev, NSW))       
     286          ALLOCATE(taulwaerstrat_mpi(klon_mpi, klev, NLW))
     287         
     288          CALL xios_send_field("tauaerstrat_in",SPREAD(tauaerstrat(:,:,:,mth_cur),1,8))
     289          CALL xios_recv_field("tauaerstrat_out",tauaerstrat_mpi)
     290          CALL xios_send_field("pizaerstrat_in",SPREAD(pizaerstrat(:,:,:,mth_cur),1,8))
     291          CALL xios_recv_field("pizaerstrat_out",pizaerstrat_mpi)
     292          CALL xios_send_field("cgaerstrat_in",SPREAD(cgaerstrat(:,:,:,mth_cur),1,8))
     293          CALL xios_recv_field("cgaerstrat_out",cgaerstrat_mpi)
     294          CALL xios_send_field("taulwaerstrat_in",SPREAD(taulwaerstrat(:,:,:,mth_cur),1,8))
     295          CALL xios_recv_field("taulwaerstrat_out",taulwaerstrat_mpi)
     296        ELSE
     297          ALLOCATE(tauaerstrat_mpi(0, 0, 0))
     298          ALLOCATE(pizaerstrat_mpi(0, 0, 0))
     299          ALLOCATE(cgaerstrat_mpi(0, 0, 0))       
     300          ALLOCATE(taulwaerstrat_mpi(0, 0, 0))
     301        ENDIF 
     302       
     303        CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
     304        CALL scatter_omp(pizaerstrat_mpi,piz_aer_strat)
     305        CALL scatter_omp(cgaerstrat_mpi,cg_aer_strat)
     306        CALL scatter_omp(taulwaerstrat_mpi,taulw_aer_strat)
     307      ELSE 
     308       
    257309!--scatter on all proc
    258       CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
    259       CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
    260       CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
    261       CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
     310        CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
     311        CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
     312        CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
     313        CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
     314
     315        DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
     316        DEALLOCATE(taulwaerstrat_mois)
     317      ENDIF
    262318
    263319      IF (is_mpi_root.AND.is_omp_root) THEN
    264 !
    265         DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat)
    266         DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
    267         DEALLOCATE(taulwaerstrat,taulwaerstrat_mois)
    268 !
    269       ENDIF !--is_mpi_root and is_omp_root
    270 
    271       DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo)
    272       DEALLOCATE(taulwaerstrat_mois_glo)
     320        DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat,taulwaerstrat)
     321      ENDIF
     322     
    273323
    274324!$OMP BARRIER

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/surf_land_mod.F90

    r3411 r3413  
    3838    USE surf_land_orchidee_nofrein_mod
    3939#else
     40#if ORCHIDEE_NOUNSTRUCT
     41    ! Compilation with cpp key ORCHIDEE_NOUNSTRUCT
     42    USE surf_land_orchidee_nounstruct_mod
     43#else
    4044    USE surf_land_orchidee_mod
     45#endif
    4146#endif
    4247#endif

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/surf_land_orchidee_mod.F90

    r3411 r3413  
    44#ifndef ORCHIDEE_NOZ0H
    55#ifndef ORCHIDEE_NOFREIN
     6#ifndef ORCHIDEE_NOUNSTRUCT
    67!
    78! This module controles the interface towards the model ORCHIDEE.
    … …  
    2324  USE cpl_mod,      ONLY : cpl_send_land_fields
    2425  USE surface_data, ONLY : type_ocean
    25   USE geometry_mod, ONLY : dx, dy
     26  USE geometry_mod, ONLY : dx, dy, boundslon, boundslat,longitude, latitude, cell_area,  ind_cell_glo
    2627  USE mod_grid_phy_lmdz
    2728  USE mod_phys_lmdz_para, mpi_root_rank=>mpi_master
    28 
     29  USE nrtype, ONLY : PI
     30 
    2931  IMPLICIT NONE
    3032
    … …  
    170172    REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: lalo
    171173    !$OMP THREADPRIVATE(lalo)
     174! boundaries of cells
     175    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: bounds_lalo
     176    !$OMP THREADPRIVATE(bounds_lalo)
    172177! pts voisins
    173178    INTEGER,ALLOCATABLE, DIMENSION(:,:), SAVE :: neighbours
    … …  
    183188    !$OMP THREADPRIVATE(lon_scat,lat_scat)
    184189
     190! area of cells
     191    REAL, ALLOCATABLE, DIMENSION (:), SAVE  :: area 
     192    !$OMP THREADPRIVATE(area)
     193
    185194    LOGICAL, SAVE                             :: lrestart_read = .TRUE.
    186195    !$OMP THREADPRIVATE(lrestart_read)
    … …  
    214223    !$OMP THREADPRIVATE(riverflow)
    215224   
     225    INTEGER :: orch_mpi_rank
     226    INTEGER :: orch_mpi_size
    216227    INTEGER :: orch_omp_rank
    217228    INTEGER :: orch_omp_size
     229
     230    REAL, ALLOCATABLE, DIMENSION(:)         :: longitude_glo
     231    REAL, ALLOCATABLE, DIMENSION(:)         :: latitude_glo
     232    REAL, ALLOCATABLE, DIMENSION(:,:)       :: boundslon_glo
     233    REAL, ALLOCATABLE, DIMENSION(:,:)       :: boundslat_glo
     234    INTEGER, ALLOCATABLE, DIMENSION(:)      :: ind_cell_glo_glo
     235    INTEGER, ALLOCATABLE, SAVE,DIMENSION(:) :: ind_cell
     236    !$OMP THREADPRIVATE(ind_cell)
     237    INTEGER :: begin, end
    218238!
    219239! Fin definition
    … …  
    258278       jg(klon) = nbp_lat
    259279
    260        IF ((.NOT. ALLOCATED(lalo))) THEN
    261           ALLOCATE(lalo(knon,2), stat = error)
     280       IF ((.NOT. ALLOCATED(area))) THEN
     281          ALLOCATE(area(knon), stat = error)
    262282          IF (error /= 0) THEN
     283             abort_message='Pb allocation area'
     284             CALL abort_physic(modname,abort_message,1)
     285          ENDIF
     286       ENDIF
     287       DO igrid = 1, knon
     288          area(igrid) = cell_area(knindex(igrid))
     289       ENDDO
     290       
     291       IF (grid_type==unstructured) THEN
     292
     293
     294         IF ((.NOT. ALLOCATED(lon_scat))) THEN
     295            ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
     296            IF (error /= 0) THEN
     297               abort_message='Pb allocation lon_scat'
     298               CALL abort_physic(modname,abort_message,1)
     299            ENDIF
     300         ENDIF
     301 
     302         IF ((.NOT. ALLOCATED(lat_scat))) THEN
     303            ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
     304            IF (error /= 0) THEN
     305               abort_message='Pb allocation lat_scat'
     306               CALL abort_physic(modname,abort_message,1)
     307            ENDIF
     308         ENDIF
     309         CALL Gather(rlon,rlon_g)
     310         CALL Gather(rlat,rlat_g)
     311
     312         IF (is_mpi_root) THEN
     313            index = 1
     314            DO jj = 2, nbp_lat-1
     315               DO ij = 1, nbp_lon
     316                  index = index + 1
     317                  lon_scat(ij,jj) = rlon_g(index)
     318                  lat_scat(ij,jj) = rlat_g(index)
     319               ENDDO
     320            ENDDO
     321            lon_scat(:,1) = lon_scat(:,2)
     322            lat_scat(:,1) = rlat_g(1)
     323            lon_scat(:,nbp_lat) = lon_scat(:,2)
     324            lat_scat(:,nbp_lat) = rlat_g(klon_glo)
     325         ENDIF
     326     
     327         CALL bcast(lon_scat)
     328         CALL bcast(lat_scat)
     329               
     330       ELSE IF (grid_type==regular_lonlat) THEN
     331
     332         IF ((.NOT. ALLOCATED(lalo))) THEN
     333            ALLOCATE(lalo(knon,2), stat = error)
     334            IF (error /= 0) THEN
     335               abort_message='Pb allocation lalo'
     336               CALL abort_physic(modname,abort_message,1)
     337            ENDIF
     338         ENDIF
     339       
     340         IF ((.NOT. ALLOCATED(bounds_lalo))) THEN
     341           ALLOCATE(bounds_lalo(knon,nvertex,2), stat = error)
     342           IF (error /= 0) THEN
    263343             abort_message='Pb allocation lalo'
    264344             CALL abort_physic(modname,abort_message,1)
    265           ENDIF
    266        ENDIF
    267        IF ((.NOT. ALLOCATED(lon_scat))) THEN
    268           ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
    269           IF (error /= 0) THEN
    270              abort_message='Pb allocation lon_scat'
    271              CALL abort_physic(modname,abort_message,1)
    272           ENDIF
    273        ENDIF
    274        IF ((.NOT. ALLOCATED(lat_scat))) THEN
    275           ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
    276           IF (error /= 0) THEN
    277              abort_message='Pb allocation lat_scat'
    278              CALL abort_physic(modname,abort_message,1)
    279           ENDIF
    280        ENDIF
    281        lon_scat = 0.
    282        lat_scat = 0.
    283        DO igrid = 1, knon
    284           index = knindex(igrid)
    285           lalo(igrid,2) = rlon(index)
    286           lalo(igrid,1) = rlat(index)
    287        ENDDO
    288 
    289        
    290        
    291        CALL Gather(rlon,rlon_g)
    292        CALL Gather(rlat,rlat_g)
    293 
    294        IF (is_mpi_root) THEN
    295           index = 1
    296           DO jj = 2, nbp_lat-1
    297              DO ij = 1, nbp_lon
    298                 index = index + 1
    299                 lon_scat(ij,jj) = rlon_g(index)
    300                 lat_scat(ij,jj) = rlat_g(index)
    301              ENDDO
    302           ENDDO
    303           lon_scat(:,1) = lon_scat(:,2)
    304           lat_scat(:,1) = rlat_g(1)
    305           lon_scat(:,nbp_lat) = lon_scat(:,2)
    306           lat_scat(:,nbp_lat) = rlat_g(klon_glo)
    307        ENDIF
     345           ENDIF
     346         ENDIF
     347       
     348         IF ((.NOT. ALLOCATED(lon_scat))) THEN
     349            ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
     350            IF (error /= 0) THEN
     351               abort_message='Pb allocation lon_scat'
     352               CALL abort_physic(modname,abort_message,1)
     353            ENDIF
     354         ENDIF
     355         IF ((.NOT. ALLOCATED(lat_scat))) THEN
     356            ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
     357            IF (error /= 0) THEN
     358               abort_message='Pb allocation lat_scat'
     359               CALL abort_physic(modname,abort_message,1)
     360            ENDIF
     361         ENDIF
     362         lon_scat = 0.
     363         lat_scat = 0.
     364         DO igrid = 1, knon
     365            index = knindex(igrid)
     366            lalo(igrid,2) = rlon(index)
     367            lalo(igrid,1) = rlat(index)
     368            bounds_lalo(igrid,:,2)=boundslon(index,:)*180./PI
     369            bounds_lalo(igrid,:,1)=boundslat(index,:)*180./PI
     370         ENDDO
     371
     372       
     373       
     374         CALL Gather(rlon,rlon_g)
     375         CALL Gather(rlat,rlat_g)
     376
     377         IF (is_mpi_root) THEN
     378            index = 1
     379            DO jj = 2, nbp_lat-1
     380               DO ij = 1, nbp_lon
     381                  index = index + 1
     382                  lon_scat(ij,jj) = rlon_g(index)
     383                  lat_scat(ij,jj) = rlat_g(index)
     384               ENDDO
     385            ENDDO
     386            lon_scat(:,1) = lon_scat(:,2)
     387            lat_scat(:,1) = rlat_g(1)
     388            lon_scat(:,nbp_lat) = lon_scat(:,2)
     389            lat_scat(:,nbp_lat) = rlat_g(klon_glo)
     390         ENDIF
    308391   
    309        CALL bcast(lon_scat)
    310        CALL bcast(lat_scat)
     392         CALL bcast(lon_scat)
     393         CALL bcast(lat_scat)
     394       
     395       ENDIF
    311396!
    312397! Allouer et initialiser le tableau des voisins et des fraction de continents
    313398!
    314        IF ( (.NOT.ALLOCATED(neighbours))) THEN
    315           ALLOCATE(neighbours(knon,8), stat = error)
    316           IF (error /= 0) THEN
    317              abort_message='Pb allocation neighbours'
    318              CALL abort_physic(modname,abort_message,1)
    319           ENDIF
    320        ENDIF
    321        neighbours = -1.
    322399       IF (( .NOT. ALLOCATED(contfrac))) THEN
    323400          ALLOCATE(contfrac(knon), stat = error)
    … …  
    334411
    335412
    336        CALL Init_neighbours(knon,neighbours,knindex,pctsrf(:,is_ter))
     413       IF (grid_type==regular_lonlat) THEN
     414 
     415         IF ( (.NOT.ALLOCATED(neighbours))) THEN
     416          ALLOCATE(neighbours(knon,8), stat = error)
     417          IF (error /= 0) THEN
     418             abort_message='Pb allocation neighbours'
     419             CALL abort_physic(modname,abort_message,1)
     420          ENDIF
     421         ENDIF
     422         neighbours = -1.
     423         CALL Init_neighbours(knon,neighbours,knindex,pctsrf(:,is_ter))
     424
     425       ELSE IF (grid_type==unstructured) THEN
     426 
     427         IF ( (.NOT.ALLOCATED(neighbours))) THEN
     428          ALLOCATE(neighbours(knon,12), stat = error)
     429          IF (error /= 0) THEN
     430             abort_message='Pb allocation neighbours'
     431             CALL abort_physic(modname,abort_message,1)
     432          ENDIF
     433         ENDIF
     434         neighbours = -1.
     435 
     436       ENDIF
     437         
    337438
    338439!
    … …  
    345446          ENDIF
    346447       ENDIF
    347        DO igrid = 1, knon
    348           ij = knindex(igrid)
    349           resolution(igrid,1) = dx(ij)
    350           resolution(igrid,2) = dy(ij)
    351        ENDDO
    352      
     448       
     449       IF (grid_type==regular_lonlat) THEN
     450         DO igrid = 1, knon
     451            ij = knindex(igrid)
     452            resolution(igrid,1) = dx(ij)
     453           resolution(igrid,2) = dy(ij)
     454         ENDDO
     455       ENDIF
     456       
    353457       ALLOCATE(coastalflow(klon), stat = error)
    354458       IF (error /= 0) THEN
    … …  
    401505    IF (debut) THEN
    402506       CALL Init_orchidee_index(knon,knindex,offset,ktindex)
    403        CALL Get_orchidee_communicator(orch_comm,orch_omp_size,orch_omp_rank)
     507       CALL Get_orchidee_communicator(orch_comm,orch_mpi_size,orch_mpi_rank, orch_omp_size,orch_omp_rank)
     508
     509       IF (grid_type==unstructured) THEN
     510         IF (knon==0) THEN
     511           begin=1
     512           end=0
     513         ELSE
     514           begin=offset+1
     515           end=offset+ktindex(knon)
     516         ENDIF
     517       
     518         IF (orch_mpi_rank==orch_mpi_size-1 .AND. orch_omp_rank==orch_omp_size-1) end=nbp_lon*nbp_lat
     519         
     520         ALLOCATE(lalo(end-begin+1,2))
     521         ALLOCATE(bounds_lalo(end-begin+1,nvertex,2))
     522         ALLOCATE(ind_cell(end-begin+1))
     523         
     524         ALLOCATE(longitude_glo(klon_glo))
     525         CALL gather(longitude,longitude_glo)
     526         CALL bcast(longitude_glo)
     527         lalo(:,2)=longitude_glo(begin:end)*180./PI
     528 
     529         ALLOCATE(latitude_glo(klon_glo))
     530         CALL gather(latitude,latitude_glo)
     531         CALL bcast(latitude_glo)
     532         lalo(:,1)=latitude_glo(begin:end)*180./PI
     533
     534         ALLOCATE(boundslon_glo(klon_glo,nvertex))
     535         CALL gather(boundslon,boundslon_glo)
     536         CALL bcast(boundslon_glo)
     537         bounds_lalo(:,:,2)=boundslon_glo(begin:end,:)*180./PI
     538 
     539         ALLOCATE(boundslat_glo(klon_glo,nvertex))
     540         CALL gather(boundslat,boundslat_glo)
     541         CALL bcast(boundslat_glo)
     542         bounds_lalo(:,:,1)=boundslat_glo(begin:end,:)*180./PI
     543         
     544         ALLOCATE(ind_cell_glo_glo(klon_glo))
     545         CALL gather(ind_cell_glo,ind_cell_glo_glo)
     546         CALL bcast(ind_cell_glo_glo)
     547         ind_cell(:)=ind_cell_glo_glo(begin:end)
     548         
     549       ENDIF
    404550       CALL Init_synchro_omp
    405551       
    406552       IF (knon > 0) THEN
    407553#ifdef CPP_VEGET
    408          CALL Init_intersurf(nbp_lon,nbp_lat,knon,ktindex,offset,orch_omp_size,orch_omp_rank,orch_comm)
     554         CALL Init_intersurf(nbp_lon,nbp_lat,knon,ktindex,offset,orch_omp_size,orch_omp_rank,orch_comm,grid=grid_type)
    409555#endif
    410556       ENDIF
    411557
    412558       
    413        IF (knon > 0) THEN
     559       IF (knon > 0) THEN 
    414560
    415561#ifdef CPP_VEGET
     562
    416563         CALL intersurf_initialize_gathered (itime+itau_phy-1, nbp_lon, nbp_lat, knon, ktindex, dtime, &
    417564               lrestart_read, lrestart_write, lalo, contfrac, neighbours, resolution, date0, &
    … …  
    421568               evap, fluxsens, fluxlat, coastalflow, riverflow, &
    422569               tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0m_new, &   
    423                lon_scat, lat_scat, q2m, t2m, z0h_new, nvm_orch)
     570               lon_scat, lat_scat, q2m(1:knon), t2m(1:knon), z0h_new(1:knon), nvm_orch, &
     571               grid=grid_type, bounds_latlon=bounds_lalo, cell_area=area, ind_cell_glo=ind_cell)
    424572#endif         
    425573       ENDIF
    … …  
    450598            evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
    451599            tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0m_new(1:knon), &
    452             lon_scat, lat_scat, q2m, t2m, z0h_new(1:knon),&
     600            lon_scat, lat_scat, q2m(1:knon), t2m(1:knon), z0h_new(1:knon),&
    453601            veget(1:knon,:),lai(1:knon,:),height(1:knon,:),&
    454602            coszang=yrmu0(1:knon))
    … …  
    525673!
    526674
    527   SUBROUTINE Get_orchidee_communicator(orch_comm,orch_omp_size,orch_omp_rank)
     675  SUBROUTINE Get_orchidee_communicator(orch_comm, orch_mpi_size, orch_mpi_rank, orch_omp_size,orch_omp_rank)
    528676  USE  mod_surf_para
    529677     
    … …  
    533681
    534682    INTEGER,INTENT(OUT) :: orch_comm
     683    INTEGER,INTENT(OUT) :: orch_mpi_size
     684    INTEGER,INTENT(OUT) :: orch_mpi_rank
    535685    INTEGER,INTENT(OUT) :: orch_omp_size
    536686    INTEGER,INTENT(OUT) :: orch_omp_rank
    … …  
    552702#ifdef CPP_MPI   
    553703      CALL MPI_COMM_SPLIT(COMM_LMDZ_PHY,color,mpi_rank,orch_comm,ierr)
     704      CALL MPI_COMM_SIZE(orch_comm,orch_mpi_size,ierr)
     705      CALL MPI_COMM_RANK(orch_comm,orch_mpi_rank,ierr)
    554706#endif
    555707   
    … …  
    683835#endif
    684836#endif
     837#endif
    685838END MODULE surf_land_orchidee_mod
Note: See TracChangeset for help on using the changeset viewer.