source: LMDZ6/trunk/libf/phylmd/readaerosolstrato.F90 @ 5087

Last change on this file since 5087 was 5084, checked in by Laurent Fairhead, 12 months ago

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, LF
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 7.8 KB
Line 
1subroutine readaerosolstrato(debut)
2
3    use netcdf95, only: nf95_close, nf95_gw_var, nf95_inq_dimid, &
4                        nf95_inq_varid, nf95_open
5    use netcdf, only: nf90_get_var, nf90_noerr, nf90_nowrite
6
7    USE phys_cal_mod, ONLY : mth_cur
8    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, &
9                                 grid2dto1d_glo, grid_type, unstructured
10    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
11    USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
12    USE mod_phys_lmdz_para
13    USE phys_state_var_mod
14    USE phys_local_var_mod
15    USE aero_mod
16    USE dimphy
17    USE print_control_mod, ONLY: prt_level,lunout
18    USE lmdz_xios
19    implicit none
20
21    include "YOMCST.h"
22
23! Variable input
24    logical debut
25
26! Variables locales
27    integer n_lat   ! number of latitudes in the input data
28    integer n_lon   ! number of longitudes in the input data
29    integer n_lev   ! number of levels in the input data
30    integer n_month ! number of months in the input data
31    real, allocatable:: latitude(:)
32    real, allocatable:: longitude(:)
33    real, allocatable:: time(:)
34    real, allocatable:: lev(:)
35    integer i, k, band, wave
36    integer, save :: mth_pre=1
37!$OMP THREADPRIVATE(mth_pre)
38
39    real, allocatable, dimension(:,:), save :: tau_aer_strat
40!$OMP THREADPRIVATE(tau_aer_strat)
41
42! Champs reconstitues
43    real, allocatable:: tauaerstrat(:, :, :, :)
44    real, allocatable:: tauaerstrat_mois(:, :, :)
45    real, allocatable:: tauaerstrat_mois_glo(:, :)
46    real, allocatable:: tau_aer_strat_mpi(:, :)
47
48! For NetCDF:
49    integer ncid_in  ! IDs for input files
50    integer varid, ncerr
51
52! Stratospheric aerosols optical properties
53! alpha_strat over the 2 bands is normalised by the 550 nm extinction coefficient
54! alpha_strat_wave is *not* normalised by the 550 nm extinction coefficient
55    real, dimension(nbands) :: alpha_strat, piz_strat, cg_strat
56    data alpha_strat/0.9922547, 0.7114912 /
57    data piz_strat  /0.9999998, 0.99762493/
58    data cg_strat   /0.73107845,0.73229635/
59    real, dimension(nwave_sw) :: alpha_strat_wave
60    data alpha_strat_wave/3.36780953,3.34667683,3.20444202,3.0293026,2.82108808/
61
62    CHARACTER (len = 20)                      :: modname = 'readaerosolstrato'
63    CHARACTER (len = 80)                      :: abort_message
64
65!--------------------------------------------------------
66
67    IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev))
68
69!--only read file if beginning of run or start of new month
70    IF (debut.OR.mth_cur.NE.mth_pre) THEN
71
72!--only root reads
73    IF (is_mpi_root.AND.is_omp_root) THEN
74
75    IF (nbands.NE.2) THEN
76        abort_message='nbands doit etre egal a 2 dans readaerosolstrat'
77        CALL abort_physic(modname,abort_message,1)
78    ENDIF
79
80    CALL nf95_open("taustrat.nc", nf90_nowrite, ncid_in)
81
82    CALL nf95_inq_varid(ncid_in, "LEV", varid)
83    CALL nf95_gw_var(ncid_in, varid, lev)
84    n_lev = size(lev)
85    IF (n_lev.NE.klev) THEN
86       abort_message='Le nombre de niveaux n est pas egal a klev'
87       CALL abort_physic(modname,abort_message,1)
88    ENDIF
89
90    CALL nf95_inq_varid(ncid_in, "LAT", varid)
91    CALL nf95_gw_var(ncid_in, varid, latitude)
92    n_lat = size(latitude)
93    WRITE(lunout,*) 'LAT aerosol strato=', n_lat, latitude
94    IF (grid_type/=unstructured) THEN
95      IF (n_lat.NE.nbp_lat) THEN
96         abort_message='Le nombre de lat n est pas egal a nbp_lat'
97         CALL abort_physic(modname,abort_message,1)
98      ENDIF
99    ENDIF
100   
101    CALL nf95_inq_varid(ncid_in, "LON", varid)
102    CALL nf95_gw_var(ncid_in, varid, longitude)
103    n_lon = size(longitude)
104    IF (grid_type/=unstructured) THEN
105      WRITE(lunout,*) 'LON aerosol strato=', n_lon, longitude
106      IF (n_lon.NE.nbp_lon) THEN
107         abort_message='Le nombre de lon n est pas egal a nbp_lon'
108         CALL abort_physic(modname,abort_message,1)
109      ENDIF
110    ENDIF
111   
112    CALL nf95_inq_varid(ncid_in, "TIME", varid)
113    CALL nf95_gw_var(ncid_in, varid, time)
114    n_month = size(time)
115    WRITE(lunout,*) 'TIME aerosol strato=', n_month, time
116    IF (n_month.NE.12) THEN
117       abort_message='Le nombre de month n est pas egal a 12'
118       CALL abort_physic(modname,abort_message,1)
119    ENDIF
120
121    IF (.not.ALLOCATED(tauaerstrat))          ALLOCATE(tauaerstrat(n_lon, n_lat, n_lev, n_month))
122    IF (.not.ALLOCATED(tauaerstrat_mois))     ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev))
123    IF (.not.ALLOCATED(tauaerstrat_mois_glo)) ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev))
124
125!--reading stratospheric AOD at 550 nm
126    CALL nf95_inq_varid(ncid_in, "TAUSTRAT", varid)
127    ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
128    WRITE(lunout,*) 'code erreur readaerosolstrato=', ncerr, varid
129
130    CALL nf95_close(ncid_in)
131
132!---select the correct month
133    IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
134     WRITE(lunout,*) 'probleme avec le mois dans readaerosolstrat =', mth_cur
135    ENDIF
136    tauaerstrat_mois(:,:,:) = tauaerstrat(:,:,:,mth_cur)
137
138!---reduce to a klon_glo grid
139    CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
140
141    ELSE
142      ALLOCATE(tauaerstrat_mois(0,0,0))
143    ENDIF !--is_mpi_root and is_omp_root
144
145!$OMP BARRIER
146
147    IF (grid_type==unstructured) THEN
148      IF (is_omp_master) THEN
149        CALL xios_send_field("taustrat_in",tauaerstrat_mois)
150        ALLOCATE(tau_aer_strat_mpi(klon_mpi, klev))
151        CALL xios_recv_field("taustrat_out",tau_aer_strat_mpi)
152      ELSE
153        ALLOCATE(tau_aer_strat_mpi(0,0))
154      ENDIF
155      CALL scatter_omp(tau_aer_strat_mpi,tau_aer_strat)
156    ELSE 
157!--scatter on all proc
158      CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
159    ENDIF
160
161!--keep memory of previous month
162    mth_pre=mth_cur
163!
164    IF (is_mpi_root.AND.is_omp_root) THEN
165!
166    DEALLOCATE(tauaerstrat)
167    DEALLOCATE(tauaerstrat_mois)
168    DEALLOCATE(tauaerstrat_mois_glo)
169!
170    ENDIF !-is_mpi_root and is_omp_root
171
172!$OMP BARRIER
173
174    ENDIF !--debut ou nouveau mois
175
176!--total vertical aod at the 6 wavelengths
177    DO wave=1, nwave_sw
178    DO k=1, klev
179    tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+tau_aer_strat(:,k)*alpha_strat_wave(wave)/alpha_strat_wave(2)
180    ENDDO
181    ENDDO
182
183!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
184    DO band=1, nbands
185!--anthropogenic aerosols bands 1 and 2
186    cg_aero(:,:,3,band)  = ( cg_aero(:,:,3,band)*piz_aero(:,:,3,band)*tau_aero(:,:,3,band) +          &
187                             cg_strat(band)*piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:) ) /  &
188                             MAX( piz_aero(:,:,3,band)*tau_aero(:,:,3,band) +                         &
189                                  piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:), 1.e-15 )
190    piz_aero(:,:,3,band)  = ( piz_aero(:,:,3,band)*tau_aero(:,:,3,band) +                             &
191                              piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:) ) /                &
192                              MAX( tau_aero(:,:,3,band) + alpha_strat(band)*tau_aer_strat(:,:), 1.e-15 )
193    tau_aero(:,:,3,band)  = tau_aero(:,:,3,band) + alpha_strat(band)*tau_aer_strat(:,:)
194!--natural aerosols bands 1 and 2
195    cg_aero(:,:,2,band)  = ( cg_aero(:,:,2,band)*piz_aero(:,:,2,band)*tau_aero(:,:,2,band) +          &
196                             cg_strat(band)*piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:) ) /  &
197                             MAX( piz_aero(:,:,2,band)*tau_aero(:,:,2,band) +                         &
198                                  piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:), 1.e-15 )
199    piz_aero(:,:,2,band)  = ( piz_aero(:,:,2,band)*tau_aero(:,:,2,band) +                             &
200                              piz_strat(band)*alpha_strat(band)*tau_aer_strat(:,:) ) /                &
201                              MAX( tau_aero(:,:,2,band) + alpha_strat(band)*tau_aer_strat(:,:),1.e-15 )
202    tau_aero(:,:,2,band)  = tau_aero(:,:,2,band) + alpha_strat(band)*tau_aer_strat(:,:)
203    ENDDO
204
205end subroutine readaerosolstrato
Note: See TracBrowser for help on using the repository browser.