Ignore:
Timestamp:
Feb 25, 2026, 1:12:05 PM (5 weeks ago)
Author:
emillour
Message:

Generic PCM:
More tidying concerning OpenMP and radiative transfert routines. Ideally all
saved variables should be threadprivate (even though it is not always necessary)
It is also better practice to only read in data from file via the master core
and then broadcast the information (works for both MPI and OpenMP).
EM

Location:
trunk/LMDZ.GENERIC/libf/phygeneric
Files:
5 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phygeneric/physiq_mod.F90

    r4081 r4083  
    2626      use rad_correlatedk_init_stellar_mod, only: rad_correlatedk_init_stellar
    2727      use rad_correlatedk_init_thermal_mod, only: rad_correlatedk_init_thermal
     28      use rad_correlatedk_read_opacity_tables_mod, only: rad_correlatedk_read_opacity_tables
    2829      use aerosol_radius, only: aerosol_radius_h2o_liquid_ice_mixture, aerosol_radius_co2
    2930      use aerosol_global_variables , only: aerosol_init, iaero_co2, iaero_h2o

  • trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_continuum_interpolation.F90

    r4077 r4083  
    2121
    2222      use datafile_mod, only: datadir
    23       use mod_phys_lmdz_para, only : is_master
     23      use mod_phys_lmdz_para, only : is_master, bcast
    2424
    2525      use gases_h, only: ngasmx, gnom, &
    … …  
    7474      double precision,save,dimension(:,:),allocatable :: abs_arr_N2N2_IR
    7575      double precision,save,dimension(:,:),allocatable :: abs_arr_N2N2_VI
    76 ! None of these saved variables are THREADPRIVATE because read by master
    77 ! and then only accessed but never modified and thus can be shared
     76!$OMP THREADPRIVATE(num_T_N2N2,temp_arr_N2N2,abs_arr_N2N2_IR,abs_arr_N2N2_VI)
    7877
    7978      ! Temperature array, continuum absorption grid for the pair O2-O2
    … …  
    8281      double precision,save,dimension(:,:),allocatable :: abs_arr_O2O2_IR
    8382      double precision,save,dimension(:,:),allocatable :: abs_arr_O2O2_VI
    84 ! None of these saved variables are THREADPRIVATE because read by master
    85 ! and then only accessed but never modified and thus can be shared
     83!$OMP THREADPRIVATE(num_T_O2O2,temp_arr_O2O2,abs_arr_O2O2_IR,abs_arr_O2O2_VI)
    8684
    8785      ! Temperature array, continuum absorption grid for the pair H2-H2
    … …  
    9088      double precision,save,dimension(:,:),allocatable :: abs_arr_H2H2_IR
    9189      double precision,save,dimension(:,:),allocatable :: abs_arr_H2H2_VI
    92 ! None of these saved variables are THREADPRIVATE because read by master
    93 ! and then only accessed but never modified and thus can be shared
     90!$OMP THREADPRIVATE(num_T_H2H2,temp_arr_H2H2,abs_arr_H2H2_IR,abs_arr_H2H2_VI)
    9491
    9592      ! Temperature array, continuum absorption grid for the pair CO2-CO2
    … …  
    9895      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CO2_IR
    9996      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CO2_VI
    100 ! None of these saved variables are THREADPRIVATE because read by master
    101 ! and then only accessed but never modified and thus can be shared
     97!$OMP THREADPRIVATE(num_T_CO2CO2,temp_arr_CO2CO2,abs_arr_CO2CO2_IR,abs_arr_CO2CO2_VI)
    10298
    10399      ! Temperature array, continuum absorption grid for the pair CH4-CH4
    … …  
    106102      double precision,save,dimension(:,:),allocatable :: abs_arr_CH4CH4_IR
    107103      double precision,save,dimension(:,:),allocatable :: abs_arr_CH4CH4_VI
    108 ! None of these saved variables are THREADPRIVATE because read by master
    109 ! and then only accessed but never modified and thus can be shared
     104!$OMP THREADPRIVATE(num_T_CH4CH4,temp_arr_CH4CH4,abs_arr_CH4CH4_IR,abs_arr_CH4CH4_VI)
    110105
    111106      ! Temperature array, continuum absorption grid for the pair H2O-H2O
    … …  
    114109      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OH2O_IR
    115110      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OH2O_VI
    116 ! None of these saved variables are THREADPRIVATE because read by master
    117 ! and then only accessed but never modified and thus can be shared
     111!$OMP THREADPRIVATE(num_T_H2OH2O,temp_arr_H2OH2O,abs_arr_H2OH2O_IR,abs_arr_H2OH2O_VI)
    118112
    119113      ! Temperature array, continuum absorption grid for the pair H2-He
    … …  
    122116      double precision,save,dimension(:,:),allocatable :: abs_arr_H2He_IR
    123117      double precision,save,dimension(:,:),allocatable :: abs_arr_H2He_VI
    124 ! None of these saved variables are THREADPRIVATE because read by master
    125 ! and then only accessed but never modified and thus can be shared
     118!$OMP THREADPRIVATE(num_T_H2He,temp_arr_H2He,abs_arr_H2He_IR,abs_arr_H2He_VI)
    126119
    127120      ! Temperature array, continuum absorption grid for the pair H2-CH4
    … …  
    130123      double precision,save,dimension(:,:),allocatable :: abs_arr_H2CH4_IR
    131124      double precision,save,dimension(:,:),allocatable :: abs_arr_H2CH4_VI
    132 ! None of these saved variables are THREADPRIVATE because read by master
    133 ! and then only accessed but never modified and thus can be shared
     125!$OMP THREADPRIVATE(num_T_H2CH4,temp_arr_H2CH4,abs_arr_H2CH4_IR,abs_arr_H2CH4_VI)
    134126
    135127      ! Temperature array, continuum absorption grid for the pair CO2-H2
    … …  
    138130      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2H2_IR
    139131      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2H2_VI
    140 ! None of these saved variables are THREADPRIVATE because read by master
    141 ! and then only accessed but never modified and thus can be shared
     132!$OMP THREADPRIVATE(num_T_CO2H2,temp_arr_CO2H2,abs_arr_CO2H2_IR,abs_arr_CO2H2_VI)
    142133
    143134      ! Temperature array, continuum absorption grid for the pair CO2-CH4
    … …  
    146137      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CH4_IR
    147138      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CH4_VI
    148 ! None of these saved variables are THREADPRIVATE because read by master
    149 ! and then only accessed but never modified and thus can be shared
     139!$OMP THREADPRIVATE(num_T_CO2CH4,temp_arr_CO2CH4,abs_arr_CO2CH4_IR,abs_arr_CO2CH4_VI)
    150140
    151141      ! Temperature array, continuum absorption grid for the pair N2-H2
    … …  
    154144      double precision,save,dimension(:,:),allocatable :: abs_arr_N2H2_IR
    155145      double precision,save,dimension(:,:),allocatable :: abs_arr_N2H2_VI
    156 ! None of these saved variables are THREADPRIVATE because read by master
    157 ! and then only accessed but never modified and thus can be shared
     146!$OMP THREADPRIVATE(num_T_N2H2,temp_arr_N2H2,abs_arr_N2H2_IR,abs_arr_N2H2_VI)
    158147
    159148      ! Temperature array, continuum absorption grid for the pair N2-CH4
    … …  
    162151      double precision,save,dimension(:,:),allocatable :: abs_arr_N2CH4_IR
    163152      double precision,save,dimension(:,:),allocatable :: abs_arr_N2CH4_VI
    164 ! None of these saved variables are THREADPRIVATE because read by master
    165 ! and then only accessed but never modified and thus can be shared
     153!$OMP THREADPRIVATE(num_T_N2CH4,temp_arr_N2CH4,abs_arr_N2CH4_IR,abs_arr_N2CH4_VI)
    166154
    167155      ! Temperature array, continuum absorption grid for the pair CO2-O2
    … …  
    170158      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2O2_IR
    171159      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2O2_VI
    172 ! None of these saved variables are THREADPRIVATE because read by master
    173 ! and then only accessed but never modified and thus can be shared
     160!$OMP THREADPRIVATE(num_T_CO2O2,temp_arr_CO2O2,abs_arr_CO2O2_IR,abs_arr_CO2O2_VI)
    174161
    175162      ! Temperature array, continuum absorption grid for the pair N2-O2
    … …  
    178165      double precision,save,dimension(:,:), allocatable :: abs_arr_N2O2_IR
    179166      double precision,save,dimension(:,:), allocatable :: abs_arr_N2O2_VI
    180 ! None of these saved variables are THREADPRIVATE because read by master
    181 ! and then only accessed but never modified and thus can be shared
     167!$OMP THREADPRIVATE(num_T_N2O2,temp_arr_N2O2,abs_arr_N2O2_IR,abs_arr_N2O2_VI)
    182168
    183169      ! Temperature array, continuum absorption grid for the pair H2O-N2
    … …  
    186172      double precision,save,dimension(:,:),allocatable :: abs_arr_H2ON2_IR
    187173      double precision,save,dimension(:,:),allocatable :: abs_arr_H2ON2_VI
    188 ! None of these saved variables are THREADPRIVATE because read by master
    189 ! and then only accessed but never modified and thus can be shared
     174!$OMP THREADPRIVATE(num_T_H2ON2,temp_arr_H2ON2,abs_arr_H2ON2_IR,abs_arr_H2ON2_VI)
    190175
    191176      ! Temperature array, continuum absorption grid for the pair H2O-O2
    … …  
    194179      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OO2_IR
    195180      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OO2_VI
    196 ! None of these saved variables are THREADPRIVATE because read by master
    197 ! and then only accessed but never modified and thus can be shared
     181!$OMP THREADPRIVATE(num_T_H2OO2,temp_arr_H2OO2,abs_arr_H2OO2_IR,abs_arr_H2OO2_VI)
    198182
    199183      ! Temperature array, continuum absorption grid for the pair H2O-CO2
    … …  
    202186      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OCO2_IR
    203187      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OCO2_VI
    204 ! None of these saved variables are THREADPRIVATE because read by master
    205 ! and then only accessed but never modified and thus can be shared
     188!$OMP THREADPRIVATE(num_T_H2OCO2,temp_arr_H2OCO2,abs_arr_H2OCO2_IR,abs_arr_H2OCO2_VI)
    206189
    207190
    208191      if(firstcall)then ! called by rad_correlatedk_read_opacity_tables  only
    209         if (is_master) print*,'----------------------------------------------------'
    210         if (is_master) print*,'Initialising continuum (rad_correlatedk_continuum_interpolation routine) from ', trim(filename)
    211 
    212 !$OMP MASTER
    213 
    214         open(unit=33, file=trim(filename), status="old", action="read",iostat=ios)
    215 
    216         if (ios.ne.0) then        ! file not found
    217           if (is_master) then
     192        if (is_master) then ! only the master needs read from files
     193          print*,'----------------------------------------------------'
     194          print*,'Initialising continuum (rad_correlatedk_continuum_interpolation routine) from ', trim(filename)
     195
     196          open(unit=33, file=trim(filename), status="old", action="read",iostat=ios)
     197
     198          if (ios.ne.0) then        ! file not found
    218199            write(*,*) 'Error from rad_correlatedk_continuum_interpolation routine'
    219200            write(*,*) 'Data file ',trim(filename),' not found.'
    … …  
    224205            write(*,*) 'Latest continuum data can be downloaded here:'
    225206            write(*,*) 'https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/continuum_data/'
     207            call abort_physic(rname,"missing input file",1)
    226208          endif
    227           call abort_physic(rname,"missing input file",1)
    228         endif
    229 
    230         ! We read the first line of the file to get the number of temperatures provided in the data file
    231         read(33, '(A)') line
    232 
    233         i = 1
    234         iT = 0
    235 
    236         do while (i .lt. len_trim(line))
    237           pos = index(line(i:), 'T=')
    238           if (pos == 0) exit
    239           i = i + pos
    240           iT = iT + 1
    241           read(line(i+2:i+10), '(E9.2)') temp_value
    242         end do
    243 
    244         num_T=iT ! num_T is the number of temperatures provided in the data file
    245        
    246         ! We read all the remaining lines of the file to get the number of wavenumbers provided in the data file
    247         iW = 0
    248         do
    249           read(33,*, end=501) line
    250           iW = iW + 1
    251         end do
     209
     210          ! We read the first line of the file to get the number of temperatures provided in the data file
     211          read(33, '(A)') line
     212
     213          i = 1
     214          iT = 0
     215
     216          do while (i .lt. len_trim(line))
     217            pos = index(line(i:), 'T=')
     218            if (pos == 0) exit
     219            i = i + pos
     220            iT = iT + 1
     221            read(line(i+2:i+10), '(E9.2)') temp_value
     222          end do
     223
     224          num_T=iT ! num_T is the number of temperatures provided in the data file
     225       
     226          ! We read all the remaining lines of the file to get the number of wavenumbers provided in the data file
     227          iW = 0
     228          do
     229            read(33,*, end=501) line
     230            iW = iW + 1
     231          end do
    252232       
    253233501 continue
    254234       
    255         num_wn=iW ! num_wn is the number of wavenumbers provided in the data file
    256        
    257         close(33)
    258 
     235          num_wn=iW ! num_wn is the number of wavenumbers provided in the data file
     236       
     237          close(33)
     238       
     239        endif ! of if (is_master)
     240       
     241        ! broadcast information to all cores
     242        call bcast(num_T)
     243        call bcast(num_wn)
     244
     245        ! allocate arrays
    259246        allocate(temp_arr(num_T))
    260247        allocate(wn_arr(num_wn))
    … …  
    263250        ! We now open and read the file a second time to extract the temperature array, wavenumber array and continuum absorption data
    264251
    265         open(unit=33, file=trim(filename), status="old", action="read")
    266        
    267         ! We extract the temperature array (temp_arr)
    268        
    269         read(33, '(A)') line
    270 
    271         i = 1
    272         iT = 0
    273 
    274         do while (i .lt. len_trim(line))
    275           pos = index(line(i:), 'T=')
    276           if (pos == 0) exit
    277           i = i + pos
    278           iT = iT + 1
    279           read(line(i+2:i+10), '(E9.2)') temp_arr(iT)
    280         end do
    281        
    282         ! We extract the wavenumber array (wn_arr) and continuum absorption (abs_arr)
    283 
    284         do iW=1,num_wn
    285           read(33,*) wn_arr(iW), (abs_arr(iW, iT), iT=1,num_T)
    286         end do
    287 
    288         close(33)
    289        
    290         print*,'We read continuum absorption data for the pair ', trim(gnom(igas_X)),'-',trim(gnom(igas_Y))
    291         print*,'Temperature grid of the dataset: ', temp_arr(:)
     252        if (is_master) then ! only the master needs to read from files
     253          open(unit=33, file=trim(filename), status="old", action="read")
     254       
     255          ! We extract the temperature array (temp_arr)
     256       
     257          read(33, '(A)') line
     258
     259          i = 1
     260          iT = 0
     261
     262          do while (i .lt. len_trim(line))
     263            pos = index(line(i:), 'T=')
     264            if (pos == 0) exit
     265            i = i + pos
     266            iT = iT + 1
     267            read(line(i+2:i+10), '(E9.2)') temp_arr(iT)
     268          end do
     269       
     270          ! We extract the wavenumber array (wn_arr) and continuum absorption (abs_arr)
     271
     272          do iW=1,num_wn
     273            read(33,*) wn_arr(iW), (abs_arr(iW, iT), iT=1,num_T)
     274          end do
     275
     276          close(33)
     277        endif ! of if (is_master)
     278
     279        ! broadcast information to all cores
     280        call bcast(temp_arr)
     281        call bcast(wn_arr)
     282        call bcast(abs_arr)
     283       
     284        if (is_master) then
     285         print*,'We read continuum absorption data for the pair ', trim(gnom(igas_X)),'-',trim(gnom(igas_Y))
     286         print*,'Temperature grid of the dataset: ', temp_arr(:)
     287        endif
    292288       
    293289        ! We loop on all molecular pairs with available continuum data and fill the corresponding array
    … …  
    458454       
    459455
    460 !$OMP END MASTER
    461 !$OMP BARRIER
    462 
    463 
    464456      endif ! firstcall
    465457

  • trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_read_opacity_tables.F90

    r4077 r4083  
     1module rad_correlatedk_read_opacity_tables_mod
     2
     3implicit none
     4
     5contains
     6
    17      subroutine rad_correlatedk_read_opacity_tables
    28
    … …  
    3945                corrk_recombin, use_premix, nrecomb_tot, rcb2tot_idx
    4046      use rad_correlatedk_continuum_interpolation_mod, only: rad_correlatedk_continuum_interpolation
     47      use util_lagrange_interpolation_mod, only: lagrange4
     48      use mod_phys_lmdz_para, only : is_master, bcast
    4149      implicit none
    4250
    … …  
    5159
    5260      ! ALLOCATABLE ARRAYS -- AS 12/2011
    53       REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE,SAVE :: gasi8, gasv8    !read by master
    54       character*20,allocatable,DIMENSION(:),SAVE :: gastype ! for check with gnom, read by master
     61      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi8, gasv8         !read by master
     62      character*20,allocatable,DIMENSION(:) :: gastype ! for check with gnom, read by master
    5563
    5664      real*8 x, xi(4), yi(4), ans, p
    … …  
    6573      if (.not. moderntracdef) use_premix=.true. ! Added by JVO for compatibility with 'old' traceur.def
    6674     
    67 !$OMP MASTER
    6875      if (use_premix) then ! use_premix flag added by JVO, thus if pure recombining then premix is skipped
    6976
     77       if (is_master) then ! only the master needs read the file
    7078!=======================================================================
    71 !     Load variable species data, exit if we have wrong database
    72       file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
    73       file_path=TRIM(datadir)//TRIM(file_id)
    74 
    75       ! check that the file exists
    76       inquire(FILE=file_path,EXIST=file_ok)
    77       if(.not.file_ok) then
     79!       Load variable species data, exit if we have wrong database
     80        file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
     81        file_path=TRIM(datadir)//TRIM(file_id)
     82
     83        ! check that the file exists
     84        inquire(FILE=file_path,EXIST=file_ok)
     85        if(.not.file_ok) then
    7886         write(*,*)'The file ',TRIM(file_path)
    7987         write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90, exiting.'
    … …  
    8391         write(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
    8492         call abort_physic("rad_correlatedk_read_opacity_tables ", "Unable to read file", 1)
    85       endif
    86 
    87       ! check that database matches varactive toggle
    88       open(111,file=TRIM(file_path),form='formatted')
    89       read(111,*) ngas
    90 
    91       if(moderntracdef) then
     93        endif
     94
     95        ! check that database matches varactive toggle
     96        open(111,file=TRIM(file_path),form='formatted')
     97        read(111,*) ngas
     98
     99        if(moderntracdef) then
    92100           ! JVO 20 - TODO : Sanity check with nspcrad !
    93101           print*, 'Warning : Sanity check on # of gases still not implemented !!'
    94       else
    95         if(ngas.ne.ngasmx)then
     102        else
     103         if(ngas.ne.ngasmx)then
    96104           print*,'Number of gases in radiative transfer data (',ngas,') does not', &
    97105                  'match that in gases.def (',ngasmx,'), exiting.'
    98106           call abort_physic("rad_correlatedk_read_opacity_tables ", "Number of gases in radiative transfer data does not match that in gases.def", 1)
    99         endif
    100       endif
    101 
    102       if(ngas.eq.1 .and. (varactive.or.varfixed))then
     107         endif
     108        endif ! of if (moderntracdef)
     109
     110        if(ngas.eq.1 .and. (varactive.or.varfixed))then
    103111         print*,'You have varactive/fixed=.true. but the database [',  &
    104112                     corrkdir(1:LEN_TRIM(corrkdir)),           &
    105113                '] has no variable species, exiting.'
    106114         call abort_physic("rad_correlatedk_read_opacity_tables ", "You have varactive/fixed=.true. but the database has no variable species",1)
    107       elseif(ngas.gt.5 .or. ngas.lt.1)then
     115        elseif(ngas.gt.5 .or. ngas.lt.1)then
    108116         print*,ngas,' species in database [',               &
    109117                     corrkdir(1:LEN_TRIM(corrkdir)),           &
    110118                '], radiative code cannot handle this.'
    111119         call abort_physic("rad_correlatedk_read_opacity_tables ", "No gas or too many gases for radiative code", 1)
    112       endif
    113 
    114       ! dynamically allocate gastype and read from Q.dat
    115       IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( ngas ) )
    116 
    117       do igas=1,ngas
     120        endif
     121
     122        ! allocate gastype and read from Q.dat
     123        ALLOCATE( gastype( ngas ) )
     124
     125        do igas=1,ngas
    118126         read(111,*) gastype(igas)
    119127         if(corrk_recombin) then
    … …  
    122130           print*,'Gas ',igas,' is ',gastype(igas)
    123131         endif
    124       enddo
    125 
    126       ! get array size, load the coefficients
    127       open(111,file=TRIM(file_path),form='formatted')
    128       read(111,*) L_REFVAR
    129       IF( .NOT. ALLOCATED( wrefvar ) ) ALLOCATE( WREFVAR(L_REFVAR) )
    130       read(111,*) wrefvar
    131       close(111)
     132        enddo ! of do igas=1,ngas
     133
     134        ! get array size, load the coefficients
     135        open(111,file=TRIM(file_path),form='formatted')
     136        read(111,*) L_REFVAR
     137        ALLOCATE( WREFVAR(L_REFVAR) )
     138        read(111,*) wrefvar
     139        close(111)
     140      endif ! of if (is_master)
     141     
     142      ! share the information with all cores
     143      call bcast(ngas)
     144      call bcast(L_REFVAR)
     145      if (.not.is_master) ALLOCATE(WREFVAR(L_REFVAR))
     146      call bcast(WREFVAR)
     147     
    132148
    133149      if(L_REFVAR.gt.1 .and. (.not.varactive) .and. (.not.varfixed))then
    … …  
    143159      else
    144160        ! Check that gastype and gnom match
    145         do igas=1,ngas
     161        if (is_master) then
     162         do igas=1,ngas
    146163           print*,'Gas ',igas,' is ',trim(gnom(igas))
    147164           if (trim(gnom(igas)).ne.trim(gastype(igas))) then
    … …  
    151168              call abort_physic("rad_correlatedk_read_opacity_tables ", "Name of a gas in radiative transfer data does not match that in gases.def",1)
    152169           endif
    153         enddo
     170         enddo
     171        endif ! of if (is_master)
    154172        print*,'Confirmed gas match in radiative transfer and gases.def!'
    155173      endif
    … …  
    165183      else
    166184        L_REFVAR = 1
    167       endif ! use_premix
     185      endif ! of if (use_premix)
    168186
    169187!=======================================================================
    170188!     Set the weighting in g-space
    171189
    172       file_id='/corrk_data/' // TRIM(corrkdir) // '/g.dat'
    173       file_path=TRIM(datadir)//TRIM(file_id)
    174 
    175       ! check that the file exists
    176       inquire(FILE=file_path,EXIST=file_ok)
    177       if(.not.file_ok) then
     190      if (is_master) then ! only the master needs read the file
     191        file_id='/corrk_data/' // TRIM(corrkdir) // '/g.dat'
     192        file_path=TRIM(datadir)//TRIM(file_id)
     193
     194        ! check that the file exists
     195        inquire(FILE=file_path,EXIST=file_ok)
     196        if(.not.file_ok) then
    178197         write(*,*)'The file ',TRIM(file_path)
    179198         write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90, exiting.'
    … …  
    183202         write(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
    184203         call abort_physic("rad_correlatedk_read_opacity_tables ", "Unable to read file", 1)
    185       endif
    186      
    187       ! check the array size is correct, load the coefficients
    188       open(111,file=TRIM(file_path),form='formatted')
    189       read(111,*) L_NGAUSS
    190       IF( .NOT. ALLOCATED( gweight ) ) ALLOCATE( GWEIGHT(L_NGAUSS) )
    191       read(111,*) gweight
    192       close(111)
     204        endif
     205     
     206        ! check the array size is correct, load the coefficients
     207        open(111,file=TRIM(file_path),form='formatted')
     208        read(111,*) L_NGAUSS
     209        ALLOCATE( GWEIGHT(L_NGAUSS) )
     210        read(111,*) gweight
     211        close(111)
    193212 
    194       ! display the values
    195       print*,'Correlated-k g-space grid:'
    196       do n=1,L_NGAUSS
     213        ! display the values
     214        print*,'Correlated-k g-space grid:'
     215        do n=1,L_NGAUSS
    197216         print*,n,'.',gweight(n)
    198       end do
    199       print*,''
     217        end do
     218        print*,''
     219      endif ! of if (is_master)
     220
     221      ! share the information with all cores
     222      call bcast(L_NGAUSS)
     223      if (.not.is_master) ALLOCATE(GWEIGHT(L_NGAUSS))
     224      call bcast(GWEIGHT)
    200225
    201226!=======================================================================
    … …  
    206231! pressure
    207232
    208       file_id='/corrk_data/' // TRIM(corrkdir) // '/p.dat'
    209       file_path=TRIM(datadir)//TRIM(file_id)
    210 
    211       ! check that the file exists
    212       inquire(FILE=file_path,EXIST=file_ok)
    213       if(.not.file_ok) then
     233      if (is_master) then ! only the master needs read the file
     234        file_id='/corrk_data/' // TRIM(corrkdir) // '/p.dat'
     235        file_path=TRIM(datadir)//TRIM(file_id)
     236
     237        ! check that the file exists
     238        inquire(FILE=file_path,EXIST=file_ok)
     239        if(.not.file_ok) then
    214240         write(*,*)'The file ',TRIM(file_path)
    215241         write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90, exiting.'
    … …  
    219245         write(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
    220246         call abort_physic("rad_correlatedk_read_opacity_tables ", "Unable to read file", 1)
    221       endif
     247        endif
    222248     
    223       ! get array size, load the coefficients
    224       open(111,file=TRIM(file_path),form='formatted')
    225       read(111,*) L_NPREF
    226       IF( .NOT. ALLOCATED( pgasref ) ) ALLOCATE( PGASREF(L_NPREF) )
    227       read(111,*) pgasref
    228       close(111)
    229       L_PINT = (L_NPREF-1)*5+1
    230       IF( .NOT. ALLOCATED( pfgasref ) ) ALLOCATE( PFGASREF(L_PINT) )
    231 
    232       ! display the values
    233       print*,'Correlated-k pressure grid (mBar):'
    234       do n=1,L_NPREF
     249        ! get array size, load the coefficients
     250        open(111,file=TRIM(file_path),form='formatted')
     251        read(111,*) L_NPREF
     252        ALLOCATE( PGASREF(L_NPREF) )
     253        read(111,*) pgasref
     254        close(111)
     255        L_PINT = (L_NPREF-1)*5+1
     256        ALLOCATE( PFGASREF(L_PINT) )
     257
     258        ! display the values
     259        print*,'Correlated-k pressure grid (mBar):'
     260        do n=1,L_NPREF
    235261         print*,n,'. 1 x 10^',pgasref(n),' mBar'
    236       end do
    237       print*,''
     262        end do
     263        print*,''
     264      endif ! of if (is_master)
     265
     266      ! share the information with all cores
     267      call bcast(L_NPREF)
     268      if (.not.is_master) ALLOCATE(PGASREF(L_NPREF))
     269      call bcast(PGASREF)
     270      call bcast(L_PINT)
     271      if (.not.is_master) ALLOCATE(PFGASREF(L_PINT)) ! values computed below
    238272
    239273      ! save the min / max matrix values
    … …  
    252286! temperature
    253287
    254       file_id='/corrk_data/' // TRIM(corrkdir) // '/T.dat'
    255       file_path=TRIM(datadir)//TRIM(file_id)
    256 
    257       ! check that the file exists
    258       inquire(FILE=file_path,EXIST=file_ok)
    259       if(.not.file_ok) then
     288      if (is_master) then ! only the master needs read the file
     289        file_id='/corrk_data/' // TRIM(corrkdir) // '/T.dat'
     290        file_path=TRIM(datadir)//TRIM(file_id)
     291
     292        ! check that the file exists
     293        inquire(FILE=file_path,EXIST=file_ok)
     294        if(.not.file_ok) then
    260295         write(*,*)'The file ',TRIM(file_path)
    261296         write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90, exiting.'
    … …  
    265300         write(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
    266301         call abort_physic("rad_correlatedk_read_opacity_tables ", "Unable to read file",1)
    267       endif
    268 
    269       ! get array size, load the coefficients
    270       open(111,file=TRIM(file_path),form='formatted')
    271       read(111,*) L_NTREF
    272       IF( .NOT. ALLOCATED( tgasref ) ) ALLOCATE( TGASREF(L_NTREF) )
    273       read(111,*) tgasref
    274       close(111)
    275 
    276       ! display the values
    277       print*,'Correlated-k temperature grid:'
    278       do n=1,L_NTREF
     302        endif
     303
     304        ! get array size, load the coefficients
     305        open(111,file=TRIM(file_path),form='formatted')
     306        read(111,*) L_NTREF
     307        ALLOCATE( TGASREF(L_NTREF) )
     308        read(111,*) tgasref
     309        close(111)
     310
     311        ! display the values
     312        print*,'Correlated-k temperature grid:'
     313        do n=1,L_NTREF
    279314         print*,n,'.',tgasref(n),' K'
    280       end do
    281 
     315        end do
     316      endif ! of if (is_master)
     317
     318      ! share the information with all cores
     319      call bcast(L_NTREF)
     320      if (.not.is_master) ALLOCATE(TGASREF(L_NTREF))
     321      call bcast(TGASREF)
     322     
    282323      ! save the min / max matrix values
    283324      tgasmin = tgasref(1)
    … …  
    285326
    286327      if(corrk_recombin) then ! even if no premix we keep a L_REFVAR=1 to store precombined at firstcall (see
    287          IF(.NOT. ALLOCATED(gasi8)) ALLOCATE(gasi8(L_NTREF,L_NPREF,L_REFVAR+nrecomb_tot,L_NSPECTI,L_NGAUSS))
    288          IF(.NOT. ALLOCATED(gasv8)) ALLOCATE(gasv8(L_NTREF,L_NPREF,L_REFVAR+nrecomb_tot,L_NSPECTV,L_NGAUSS))
     328         ALLOCATE(gasi8(L_NTREF,L_NPREF,L_REFVAR+nrecomb_tot,L_NSPECTI,L_NGAUSS))
     329         ALLOCATE(gasv8(L_NTREF,L_NPREF,L_REFVAR+nrecomb_tot,L_NSPECTV,L_NGAUSS))
    289330      else
    290          IF(.NOT. ALLOCATED(gasi8)) ALLOCATE(gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS))
    291          IF(.NOT. ALLOCATED(gasv8)) ALLOCATE(gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS))
     331         ALLOCATE(gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS))
     332         ALLOCATE(gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS))
    292333      endif
    293 !$OMP END MASTER
    294 !$OMP BARRIER
    295334
    296335! JVO 20 : In these gasi/gasi8 matrix  we stack in same dim. variable specie and species to recombine (to keep code small)
    … …  
    299338! allocate the multidimensional arrays in radcommon_h
    300339     if(corrk_recombin) then
    301         IF(.NOT. ALLOCATED(gasi)) ALLOCATE(gasi(L_NTREF,L_PINT,L_REFVAR+nrecomb_tot,L_NSPECTI,L_NGAUSS))
    302         IF(.NOT. ALLOCATED(gasv)) ALLOCATE(gasv(L_NTREF,L_PINT,L_REFVAR+nrecomb_tot,L_NSPECTV,L_NGAUSS))
     340        ALLOCATE(gasi(L_NTREF,L_PINT,L_REFVAR+nrecomb_tot,L_NSPECTI,L_NGAUSS))
     341        ALLOCATE(gasv(L_NTREF,L_PINT,L_REFVAR+nrecomb_tot,L_NSPECTV,L_NGAUSS))
    303342        ! display the values
    304343        print*,''
    … …  
    306345        print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR+nrecomb_tot,' (',L_REFVAR,'+',nrecomb_tot,'),',L_NGAUSS,']'
    307346      else
    308         IF(.NOT. ALLOCATED(gasi)) ALLOCATE(gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS))
    309         IF(.NOT. ALLOCATED(gasv)) ALLOCATE(gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS))
     347        ALLOCATE(gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS))
     348        ALLOCATE(gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS))
    310349        ! display the values
    311350        print*,''
    … …  
    359398      End if
    360399
    361 !$OMP MASTER         
    362 
    363       ! VISIBLE
    364       if (callgasvis) then
     400      if (is_master) then ! only the master needs to fill gasv8
     401
     402       ! VISIBLE
     403       if (callgasvis) then
    365404           
    366405        ! Looking for premixed corrk files corrk_gcm_VI.dat if needed
    … …  
    416455        if(nVI_limit.eq.0) then
    417456         gasv8(:,:,:,:,:)= gasv8(:,:,:,:,:)+kappa_VI
    418            else if (nVI_limit.eq.L_NSPECTV) then
     457        else if (nVI_limit.eq.L_NSPECTV) then
    419458         gasv8(:,:,:,:,:)= gasv8(:,:,:,:,:)+kappa_IR
    420       else
     459        else
    421460         gasv8(:,:,:,1:nVI_limit,:)= gasv8(:,:,:,1:nVI_limit,:)+kappa_IR
    422461         gasv8(:,:,:,nVI_limit+1:L_NSPECTV,:)= gasv8(:,:,:,nVI_limit+1:L_NSPECTV,:)+kappa_VI
    423       end if
     462        end if
    424463           
    425          else
     464       else
    426465            print*,'Visible corrk gaseous absorption is set to zero.'
    427466            gasv8(:,:,:,:,:)=0.0
    428          endif ! callgasvis
    429          
    430 !$OMP END MASTER
    431 !$OMP BARRIER
    432 
     467       endif ! callgasvis
     468     
     469     endif ! of (is_master)
     470     
     471     ! share the information with all cores
     472     call bcast(gasv8)
     473     
     474      if (is_master) then ! only the master needs to fill gasi8
    433475      ! INFRA-RED
    434476     
    … …  
    437479        if ((corrkdir(1:4).eq.'null'))then       !.or.(TRIM(corrkdir).eq.'null_LowTeffStar')) then
    438480           print*,'Infrared corrk gaseous absorption is set to zero if graybody=F'
    439 !$OMP MASTER         
    440481           gasi8(:,:,:,:,:)=0.0
    441 !$OMP END MASTER
    442 !$OMP BARRIER
    443482        else
    444483           file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
    … …  
    454493           endif
    455494         
    456 !$OMP MASTER                   
    457495           open(111,file=TRIM(file_path),form='formatted')
    458496           read(111,*) gasi8(:,:,:L_REFVAR,:,:)
    459497           close(111)
    460 !$OMP END MASTER
    461 !$OMP BARRIER
    462      
     498
    463499           ! 'fzero' is a currently unused feature that allows optimisation
    464500           ! of the radiative transfer by neglecting bands where absorption
    … …  
    520556             call abort_physic("rad_correlatedk_read_opacity_tables ", "No absorption data found", 1)
    521557          endif
    522 !$OMP MASTER           
     558
    523559          open(111,file=TRIM(file_path),form='formatted')
    524560          read(111,*) gasi8(:,:,L_REFVAR+igas,:,:)
    525561          close(111)
    526 !$OMP END MASTER
    527 !$OMP BARRIER
    528562        enddo                           
    529563      endif ! corrk_recombin
    530564
    531 !$OMP MASTER                 
     565      endif ! of if (is_master)
     566
     567      ! share the information with all cores
     568      call bcast(gasi8)
     569      call bcast(fzeroi)
     570      call bcast(fzerov)
     571     
    532572      if(nIR_limit.eq.0) then
    533573         gasi8(:,:,:,:,:)= gasi8(:,:,:,:,:)+kappa_VI
    … …  
    568608         end do
    569609      end do
    570 !$OMP END MASTER
    571 !$OMP BARRIER
    572610
    573611!     Interpolate the values:  first the longwave
    … …  
    591629                     yi(3) = gasi8(nt,n+2,nh,nw,ng)
    592630                     yi(4) = gasi8(nt,n+3,nh,nw,ng)
    593                      call lagrange(x,xi,yi,ans)
     631                     call lagrange4(x,xi,yi,ans)
    594632                     gasi(nt,m,nh,nw,ng) = 10.0**ans
    595633                  end do
    … …  
    607645                        yi(3) = gasi8(nt,n+1,nh,nw,ng)
    608646                        yi(4) = gasi8(nt,n+2,nh,nw,ng)
    609                         call lagrange(x,xi,yi,ans)
     647                        call lagrange4(x,xi,yi,ans)
    610648                        gasi(nt,i,nh,nw,ng) = 10.0**ans
    611649                     end do
    … …  
    626664                     yi(3) = gasi8(nt,n,nh,nw,ng)
    627665                     yi(4) = gasi8(nt,n+1,nh,nw,ng)
    628                      call lagrange(x,xi,yi,ans)
     666                     call lagrange4(x,xi,yi,ans)
    629667                     gasi(nt,i,nh,nw,ng) = 10.0**ans
    630668                  end do 
    … …  
    660698                     yi(3) = gasv8(nt,n+2,nh,nw,ng)
    661699                     yi(4) = gasv8(nt,n+3,nh,nw,ng)
    662                      call lagrange(x,xi,yi,ans)
     700                     call lagrange4(x,xi,yi,ans)
    663701                     gasv(nt,m,nh,nw,ng) = 10.0**ans
    664702                  end do
    … …  
    676714                        yi(3) = gasv8(nt,n+1,nh,nw,ng)
    677715                        yi(4) = gasv8(nt,n+2,nh,nw,ng)
    678                         call lagrange(x,xi,yi,ans)
     716                        call lagrange4(x,xi,yi,ans)
    679717                        gasv(nt,i,nh,nw,ng) = 10.0**ans
    680718                     end do
    … …  
    695733                     yi(3) = gasv8(nt,n,nh,nw,ng)
    696734                     yi(4) = gasv8(nt,n+1,nh,nw,ng)
    697                      call lagrange(x,xi,yi,ans)
     735                     call lagrange4(x,xi,yi,ans)
    698736                     gasv(nt,i,nh,nw,ng) = 10.0**ans
    699737                  end do 
    … …  
    813851
    814852!     Deallocate local arrays
    815 !$OMP BARRIER
    816 !$OMP MASTER
    817       IF( ALLOCATED( gasi8 ) ) DEALLOCATE( gasi8 )
    818       IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
    819       IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
    820       IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
    821 !$OMP END MASTER
    822 !$OMP BARRIER
     853      DEALLOCATE( gasi8 )
     854      DEALLOCATE( gasv8 )
     855      if (is_master) DEALLOCATE( gastype )
    823856
    824857    end subroutine rad_correlatedk_read_opacity_tables
     858
     859end module rad_correlatedk_read_opacity_tables_mod

  • trunk/LMDZ.GENERIC/libf/phygeneric/radcommon_h.F90

    r4081 r4083  
    6363!
    6464
    65       REAL*8 BWNI(L_NSPECTI+1) !BWNI read by master in rad_correlatedk_init_thermal
     65      REAL*8,SAVE :: BWNI(L_NSPECTI+1) !BWNI read in rad_correlatedk_init_thermal
    6666!$OMP THREADPRIVATE(BWNI)
    67       REAL*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI)
     67      REAL*8,SAVE :: WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI)
    6868!$OMP THREADPRIVATE(WNOI,DWNI,WAVEI)
    69       REAL*8 BWNV(L_NSPECTV+1) !BWNV read by master in rad_correlatedk_init_stellar
     69      REAL*8,SAVE :: BWNV(L_NSPECTV+1) !BWNV read by master in rad_correlatedk_init_stellar
    7070!$OMP THREADPRIVATE(BWNV)
    71       REAL*8 WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV)
    72       REAL*8 STELLARF(L_NSPECTV)
     71      REAL*8,SAVE :: WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV)
     72      REAL*8,SAVE :: STELLARF(L_NSPECTV)
    7373!$OMP THREADPRIVATE(WNOV,DWNV,WAVEV,STELLARF)
    7474
    75       REAL*8 blami(L_NSPECTI+1)
    76       REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90
     75      REAL*8,SAVE :: blami(L_NSPECTI+1)
     76      REAL*8,SAVE :: blamv(L_NSPECTV+1) ! these are needed by suaer.F90
    7777!$OMP THREADPRIVATE(blami,blamv)
    7878
    79       !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 
    80       REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv
    81       REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, WREFVAR, PFGASREF, GWEIGHT
    82       real*8 FZEROI(L_NSPECTI)
    83       real*8 FZEROV(L_NSPECTV)
    84       real*8 pgasmin, pgasmax
    85       real*8 tgasmin, tgasmax
    86 !$OMP THREADPRIVATE(gasi,gasv,&  !wrefvar,pgasref,tgasref,pfgasref read by master in rad_correlatedk_read_opacity_tables
    87         !$OMP FZEROI,FZEROV)     !pgasmin,pgasmax,tgasmin,tgasmax read by master in rad_correlatedk_read_opacity_tables
     79!gasi and gasv initialized in rad_correlatedk_read_opacity_tables
     80      REAL*8,DIMENSION(:,:,:,:,:),ALLOCATABLE,SAVE :: gasi, gasv
     81!$OMP THREADPRIVATE(gasi,gasv)
     82!PGASREF, TGASREF, WREFVAR initialized in rad_correlatedk_read_opacity_tables
     83      REAL*8,DIMENSION(:),ALLOCATABLE,SAVE :: PGASREF, TGASREF, WREFVAR
     84!$OMP THREADPRIVATE(PGASREF, TGASREF, WREFVAR)
     85!PFGASREF, GWEIGHT initialized in rad_correlatedk_read_opacity_tables
     86      REAL*8,DIMENSION(:),ALLOCATABLE,SAVE :: PFGASREF, GWEIGHT
     87!$OMP THREADPRIVATE(PFGASREF, GWEIGHT)
     88!FREROI and FZEROV initialized in rad_correlatedk_read_opacity_tables
     89      real*8,save :: FZEROI(L_NSPECTI)
     90      real*8,save :: FZEROV(L_NSPECTV)
     91!$OMP THREADPRIVATE(FZEROI,FZEROV)
     92!pgasmin,pgasmax,tgasmin,tgasmax read in rad_correlatedk_read_opacity_tables
     93      real*8,save :: pgasmin, pgasmax
     94      real*8,save :: tgasmin, tgasmax
     95!$OMP THREADPRIVATE(pgasmin,pgasmax,tgasmin,tgasmax)
    8896
    8997      real,allocatable,save :: QVISsQREF(:,:,:)
    … …  
    121129!      REAL :: omegaREFvis(naerkind,nsizemax)
    122130      real,allocatable,save :: omegaREFir(:,:)
     131!$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFir)
    123132
    124133      REAL,SAVE :: tstellar ! Stellar brightness temperature (SW)
    … …  
    127136
    128137      real*8,save :: PTOP
     138!$OMP THREADPRIVATE(tstellar,planckir,PTOP)
    129139
    130140      real*8,parameter :: UBARI = 0.5D0
    131141
    132 !$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFir,&         ! gweight read by master in rad_correlatedk_read_opacity_tables
    133                 !$OMP tstellar,planckir,PTOP)
    134142
    135143!     If the gas optical depth (top to the surface) is less than

  • trunk/LMDZ.GENERIC/libf/phygeneric/radinc_h.F90

    r4077 r4083  
    1010
    1111!======================================================================
    12 !
    13 !     RADINC.H
    1412!
    1513!     Includes for the radiation code; RADIATION LAYERS, LEVELS,
    … …  
    5351!               can in princple be anything: currently it's H2O.
    5452!
    55 !     NAERKIND  The number of radiatively active aerosol types
    56 !
    5753!     NSIZEMAX  The maximum number of aerosol particle sizes
    5854!
    … …  
    6460!$OMP THREADPRIVATE(L_NLAYRAD,L_LEVELS,L_NLEVRAD)
    6561
    66       ! These are set in rad_correlatedk_read_opacity_tables
    67       ! [uses allocatable arrays] -- AS 12/2011
    68       integer :: L_NPREF, L_NTREF, L_REFVAR, L_PINT, L_NGAUSS  !L_NPREF, L_NTREF, L_REFVAR, L_PINT, L_NGAUSS read by master in rad_correlatedk_read_opacity_tables
     62      ! L_NPREF, L_NTREF, L_REFVAR, L_PINT, L_NGAUSS initialized
     63      ! in rad_correlatedk_read_opacity_tables
     64      integer,save :: L_NPREF, L_NTREF, L_REFVAR, L_PINT, L_NGAUSS 
     65!$OMP THREADPRIVATE(L_NPREF,L_NTREF,L_REFVAR,L_PINT,L_NGAUSS)
    6966
    7067      integer, parameter :: L_NSPECTI = NBinfrared
    7168      integer, parameter :: L_NSPECTV = NBvisible
    7269
    73 !      integer, parameter :: NAERKIND  = 2 ! set in scatterers.h
    7470      real,    parameter :: L_TAUMAX  = 35
    7571
    … …  
    8278      !             Default is wide range : 30K-1500K, with 0.1K step
    8379      !             ->  NTstart=300, Nstop=15000, NTfac=10
    84       integer :: NTstart, NTstop
    85       real*8  :: NTfac
     80      ! NTstart, NTstop, NTfac initialized by ini_radinc_h
     81      integer,save :: NTstart, NTstop
     82      real*8,save  :: NTfac
     83!$OMP THREADPRIVATE(NTstart,NTstop,NTfac)
    8684
    8785      ! Maximum number of grain size classes for aerosol convolution:

  • trunk/LMDZ.GENERIC/libf/phygeneric/util_lagrange_interpolation.F90

    r4082 r4083  
    1       subroutine lagrange(x, xi, yi, ans)
     1module util_lagrange_interpolation_mod
    22
    3 C  Lagrange interpolation - Polynomial interpolation at point x
    4 C  xi(1) <= x <= xi(4).  Yi(n) is the functional value at XI(n).
     3implicit none
     4
     5contains
     6
     7!======================================================================!
     8      subroutine lagrange4(x, xi, yi, ans)
     9
     10!  Lagrange interpolation - Polynomial interpolation at point x
     11!  (3rd order since relying on 4 points)
     12!  xi(1) <= x <= xi(4).  Yi(n) is the functional value at XI(n).
    513
    614      implicit none
    715
    8       real*8 x, xi(4), yi(4), ans
     16      real*8,intent(in) :: x, xi(4), yi(4)
     17      real*8, intent(out) :: ans
     18
    919      real*8 fm1, fm2, fm3, fm4
    1020
    11 C======================================================================!
     21!======================================================================!
    1222
    1323      fm1   = x - XI(1)
    … …  
    1626      fm4   = x - XI(4)
    1727
    18 C  Get the answer at the requested X
     28!  Get the answer at the requested X
    1929 
    20       ans = fm2*fm3*fm4*YI(1)/
    21      *                ((XI(1)-XI(2))*(XI(1)-XI(3))*(XI(1)-XI(4)))  +
    22      *      fm1*fm3*fm4*YI(2)/
    23      *                ((XI(2)-XI(1))*(XI(2)-XI(3))*(XI(2)-XI(4)))  +
    24      *      fm1*fm2*fm4*YI(3)/
    25      *                ((XI(3)-XI(1))*(XI(3)-XI(2))*(XI(3)-XI(4)))  +
    26      *      fm1*fm2*fm3*YI(4)/
    27      *                ((XI(4)-XI(1))*(XI(4)-XI(2))*(XI(4)-XI(3)))
     30      ans = fm2*fm3*fm4*YI(1)/                                        &
     31                      ((XI(1)-XI(2))*(XI(1)-XI(3))*(XI(1)-XI(4)))  +  &
     32            fm1*fm3*fm4*YI(2)/                                        &
     33                      ((XI(2)-XI(1))*(XI(2)-XI(3))*(XI(2)-XI(4)))  +  &
     34            fm1*fm2*fm4*YI(3)/                                        &
     35                      ((XI(3)-XI(1))*(XI(3)-XI(2))*(XI(3)-XI(4)))  +  &
     36            fm1*fm2*fm3*YI(4)/                                        &
     37                      ((XI(4)-XI(1))*(XI(4)-XI(2))*(XI(4)-XI(3)))
    2838
    29       return
    30       end
     39      end subroutine lagrange4
     40
     41!======================================================================!
     42
     43end module util_lagrange_interpolation_mod
Note: See TracChangeset for help on using the changeset viewer.