Changeset 2820 for LMDZ5/trunk
- Timestamp:
- Mar 10, 2017, 5:46:08 PM (8 years ago)
- Location:
- LMDZ5/trunk/libf/phylmd
- Files:
-
- 4 edited
-
conf_phys_m.F90 (modified) (1 diff)
-
open_climoz_m.F90 (modified) (1 diff)
-
physiq_mod.F90 (modified) (2 diffs)
-
regr_pr_time_av_m.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phylmd/conf_phys_m.F90
r2788 r2820 2048 2048 ! 2049 2049 ! 2050 ok_daily_climoz = . TRUE.2050 ok_daily_climoz = .FALSE. 2051 2051 CALL getin('ok_daily_climoz', ok_daily_climoz_omp) 2052 2052 ! -
LMDZ5/trunk/libf/phylmd/open_climoz_m.F90
r2788 r2820 1 1 ! $Id$ 2 moduleopen_climoz_m2 MODULE open_climoz_m 3 3 4 implicit none 4 USE print_control_mod, ONLY: lunout 5 IMPLICIT NONE 5 6 6 contains 7 CONTAINS 7 8 8 subroutine open_climoz(ncid, press_in_cen, press_in_edg, time_in) 9 !------------------------------------------------------------------------------- 10 ! 11 SUBROUTINE open_climoz(ncid, press_in_cen, press_in_edg, time_in, daily, adjust) 12 ! 13 !------------------------------------------------------------------------------- 14 USE netcdf95, ONLY: nf95_open, nf95_close, nf95_gw_var, nf95_inq_varid 15 USE netcdf, ONLY: nf90_nowrite 16 USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root 17 USE mod_phys_lmdz_mpi_transfert, ONLY: bcast_mpi 18 USE phys_cal_mod, ONLY: calend, year_len, year_cur 19 !------------------------------------------------------------------------------- 20 ! Purpose: This procedure should be called once per "gcm" run, by a single 21 ! thread of each MPI process. 22 ! The root MPI process opens "climoz_LMDZ.nc", reads the pressure 23 ! levels and the times and broadcasts them to the other processes. 24 ! We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa 25 ! and in strictly ascending order. 26 !------------------------------------------------------------------------------- 27 ! Arguments (OpenMP shared): 28 INTEGER, INTENT(OUT):: ncid !--- "climoz_LMDZ.nc" identifier 29 REAL, POINTER :: press_in_cen(:) !--- at cells centers 30 REAL, POINTER :: press_in_edg(:) !--- at the interfaces (pressure intervals) 31 REAL, POINTER :: time_in(:) !--- records times, in days since Jan. 1st 32 LOGICAL, INTENT(IN) :: daily !--- daily files (calendar dependent days nb) 33 LOGICAL, INTENT(IN) :: adjust !--- tropopause adjustement required 34 ! pressure levels press_in_cen/edg are in Pa a,d strictly ascending order. 35 ! time_in is only used for monthly files (14 records) 36 !------------------------------------------------------------------------------- 37 ! Local variables: 38 INTEGER :: varid !--- NetCDF variables identifier 39 INTEGER :: nlev, ntim !--- pressure levels and time records numbers 40 CHARACTER(LEN=80) :: sub 41 CHARACTER(LEN=320) :: msg 42 !------------------------------------------------------------------------------- 43 sub="open_climoz" 44 WRITE(lunout,*)"Entering routine "//TRIM(sub) 9 45 10 ! This procedure should be called once per "gcm" run, by a single 11 ! thread of each MPI process. 12 ! The root MPI process opens "climoz_LMDZ.nc", reads the pressure 13 ! levels and the times and broadcasts them to the other processes. 46 IF(is_mpi_root) THEN 14 47 15 ! We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa 16 ! and in strictly ascending order. 48 !--- OPEN FILE, READ PRESSURE LEVELS AND TIME VECTOR 49 CALL nf95_open("climoz_LMDZ.nc", nf90_nowrite, ncid) 50 CALL nf95_inq_varid(ncid, "plev", varid) 51 CALL nf95_gw_var(ncid, varid, press_in_cen) 52 ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa: 53 press_in_cen = press_in_cen * 100. 54 nlev = SIZE(press_in_cen) 55 CALL NF95_INQ_VARID(ncID, "time", varID) 56 CALL NF95_GW_VAR(ncid, varid, time_in) 57 ntim = SIZE(time_in) 17 58 18 use netcdf95, only: nf95_open, nf95_close, nf95_gw_var, nf95_inq_varid 19 use netcdf, only: nf90_nowrite 59 !--- BUILD EDGES OF PRESSURE INTERVALS: HALFWAY IN LOGARITHMS 60 ALLOCATE(press_in_edg(nlev+1)) 61 press_in_edg=[0.,press_in_cen(1:nlev-1)*press_in_cen(2:nlev),HUGE(0.)] 20 62 21 use mod_phys_lmdz_mpi_data, only: is_mpi_root 22 use mod_phys_lmdz_mpi_transfert, only: bcast_mpi ! broadcast 23 USE phys_cal_mod, ONLY: calend, year_len, year_cur 63 !--- CHECK RECORDS NUMBER AND DISPLAY CORRESPONDING INFORMATION 64 IF(daily.AND.ntim/=year_len) THEN 65 WRITE(msg,'(a,3(i4,a))')TRIM(sub)//': Expecting a daily ozone file with',& 66 &year_len,' records (year ',year_cur,') ; found ',ntim,' instead' 67 CALL abort_physic(sub, msg, 1) 68 ELSE IF(ALL([360,14]/=ntim)) THEN 69 WRITE(msg,'(a,i4,a)')TRIM(sub)//': Expecting an ozone file with 14 (mont'& 70 &//'hly case) or 360 (old style files) records ; found ',ntim,' instead' 71 CALL abort_physic(sub, msg, 1) 72 ELSE 73 IF(daily) THEN 74 WRITE(msg,'(a,2(i4,a))')'daily file (',ntim,' days in ',year_cur,')' 75 ELSE IF(ntim==14) THEN 76 msg='14 records monthly file' 77 ELSE 78 msg='360 records files (old convention)' 79 END IF 80 WRITE(lunout,*)TRIM(sub)//': Using a '//TRIM(msg) 81 END IF 24 82 25 ! OpenMP shares arguments: 26 integer, intent(out):: ncid ! of "climoz_LMDZ.nc", OpenMP shared 83 !--- MESSAGE ABOUT OPTIONAL STRETCHING FOR TROPOPAUSES MATCHING 84 IF(adjust) THEN 85 WRITE(lunout,*)TRIM(sub)//': Adjusting O3 field to match gcm tropopause.' 86 ELSE 87 WRITE(lunout,*)TRIM(sub)//': Interpolating O3 field directly on gcm levels.' 88 END IF 27 89 28 ! pressure levels for ozone climatology, in Pa, strictly ascending order 29 real, pointer:: press_in_cen(:) !--- at cells centers 30 real, pointer:: press_in_edg(:) !--- at the interfaces (pressure intervals) 90 END IF 91 CALL bcast_mpi(nlev) 92 IF(.NOT.is_mpi_root) ALLOCATE(press_in_cen(nlev )); CALL bcast_mpi(press_in_cen) 93 IF(.NOT.is_mpi_root) ALLOCATE(press_in_edg(nlev+1)); CALL bcast_mpi(press_in_edg) 94 CALL bcast_mpi(ntim) 95 IF(.NOT.is_mpi_root) ALLOCATE(time_in(ntim)); CALL bcast_mpi(time_in) 31 96 32 real, pointer:: time_in(:) 33 ! times of ozone climatology records, in days since Jan. 1st 0h 97 END SUBROUTINE open_climoz 98 ! 99 !------------------------------------------------------------------------------- 34 100 35 ! Variables local to the procedure: 36 37 integer varid ! for NetCDF 38 integer n_plev ! number of pressure levels in the input data 39 integer n_time ! number of time records in the input field 40 integer k 41 CHARACTER(LEN=80) :: modname 42 CHARACTER(LEN=320) :: msg 43 44 !--------------------------------------- 45 46 modname="open_climoz" 47 print *, "Call sequence information: "//TRIM(modname) 48 49 if (is_mpi_root) then 50 call nf95_open("climoz_LMDZ.nc", nf90_nowrite, ncid) 51 52 call nf95_inq_varid(ncid, "plev", varid) 53 call nf95_gw_var(ncid, varid, press_in_cen) 54 ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa: 55 press_in_cen = press_in_cen * 100. 56 n_plev = size(press_in_cen) 57 CALL NF95_INQ_VARID(ncID, "time", varID) 58 CALL NF95_GW_VAR(ncid, varid, time_in) 59 n_time = SIZE(time_in) 60 IF(ALL([year_len,14]/=n_time)) THEN !--- Check records number 61 WRITE(msg,'(a,3(i4,a))')TRIM(modname)//': expected input file could be& 62 & monthly or daily (14 or ',year_len,' records for year ',year_cur,' a& 63 &nd calendar "'//TRIM(calend)//'"), but ',n_time,' records were found' 64 CALL abort_physic(modname, msg, 1) 65 END IF 66 67 end if 68 69 call bcast_mpi(n_plev) 70 if (.not. is_mpi_root) allocate(press_in_cen(n_plev)) 71 call bcast_mpi(press_in_cen) 72 call bcast_mpi(n_time) 73 if (.not. is_mpi_root) allocate(time_in(n_time)) 74 call bcast_mpi(time_in) 75 76 ! Compute edges of pressure intervals: 77 allocate(press_in_edg(n_plev + 1)) 78 if (is_mpi_root) then 79 press_in_edg(1) = 0. 80 ! We choose edges halfway in logarithm: 81 forall (k = 2:n_plev) press_in_edg(k) = sqrt(press_in_cen(k - 1) * press_in_cen(k)) 82 press_in_edg(n_plev + 1) = huge(0.) 83 ! (infinity, but any value guaranteed to be greater than the 84 ! surface pressure would do) 85 end if 86 call bcast_mpi(press_in_edg) 87 ! deallocate(press_in_cen) ! pointer 88 89 end subroutine open_climoz 90 91 end module open_climoz_m 101 END MODULE open_climoz_m -
LMDZ5/trunk/libf/phylmd/physiq_mod.F90
r2819 r2820 1670 1670 !$omp single 1671 1671 IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, & 1672 press_edg_climoz, time_climoz)1672 press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause) 1673 1673 !$omp end single 1674 1674 ! … … 1975 1975 1976 1976 wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz) 1977 ELSE IF(ANY([1,2]==read_climoz)) THEN1978 ro3i = REAL(INT((days_elapsed+jh_cur-jh_1jan)/year_len*360.))1979 IF(INT(ro3i) == 360) ro3i = 359.1980 ! (This should never occur, except perhaps because of roundup1981 ! error. See documentation.)1982 WRITE(*,*)' >> Using 360 records ozone files (old method).'1983 CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &1984 ro3i+1., press_edg_climoz, paprs, wo)1985 FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &1986 * zmasse / dobson_u / 1e31987 1977 ELSE 1978 !--- ro3i = elapsed days number since current year 1st january, 0h 1979 ro3i=days_elapsed+jh_cur-jh_1jan 1980 !--- scaling for old style files (360 records) 1981 IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len 1988 1982 IF(adjust_tropopause) THEN 1989 WRITE(*,*)' >> Adjusting O3 field to match gcm tropopause.' 1990 CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot, & 1991 pr_tropopause) 1992 CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz-2), & 1993 days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo, & 1994 time_climoz, latitude_deg, press_cen_climoz, pr_tropopause) 1983 CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot, & 1984 pr_tropopause) 1985 CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), & 1986 ro3i, press_edg_climoz, paprs, wo, time_climoz, & 1987 latitude_deg, press_cen_climoz, pr_tropopause) 1995 1988 ELSE 1996 WRITE(*,*)' >> Interpolating O3 field directly on gcm levels.' 1997 CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz-2), & 1998 days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo, & 1999 time_climoz) 1989 CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), & 1990 ro3i, press_edg_climoz, paprs, wo, time_climoz) 2000 1991 END IF 2001 1992 ! Convert from mole fraction of ozone to column density of ozone in a -
LMDZ5/trunk/libf/phylmd/regr_pr_time_av_m.F90
r2819 r2820 2 2 MODULE regr_pr_time_av_m 3 3 4 USE print_control_mod, ONLY: lunout 4 5 USE write_field_phy 5 6 USE mod_phys_lmdz_transfert_para, ONLY: bcast … … 117 118 CHARACTER(LEN=80) :: sub 118 119 CHARACTER(LEN=320) :: msg 119 INTEGER :: vID, ncerr, n_var, nlev_in, ndim, i, ibot, itop, itrp,itrp0120 LOGICAL :: l trop!--- Need to adjust tropopause120 INTEGER :: vID, ncerr, n_var, nlev_in,ntim_in, ndim, i, ibot, itop, itrp,itrp0 121 LOGICAL :: lAdjTro !--- Need to adjust tropopause 121 122 REAL :: y_frac !--- Elapsed year fraction 122 123 REAL :: alpha, beta, al !--- For strectching/interpolation … … 137 138 !------------------------------------------------------------------------------- 138 139 sub="regr_pr_time_av" 139 nlev_in=SIZE(pint_in)-1 140 nlev_in=SIZE(pint_in)-1; ntim_in=SIZE(time_in) 140 141 CALL assert(SIZE(v3,1)==klon, TRIM(sub)//" v3 klon") 141 142 CALL assert(SIZE(v3,2)==nbp_lev, TRIM(sub)//" v3 nbp_lev") … … 145 146 IF(PRESENT(ptrop_ou)) CALL assert(SIZE(ptrop_ou)==klon,TRIM(sub)//" ptrop_ou klon") 146 147 IF(PRESENT(pcen_in)) CALL assert(SIZE(pcen_in)==nlev_in,TRIM(sub)//" pcen_in") 147 ltrop=PRESENT(lat_in).AND.PRESENT(pcen_in).AND.PRESENT(ptrop_ou) 148 lAdjTro=PRESENT(ptrop_ou) 149 IF(lAdjTro.AND.(.NOT.PRESENT(lat_in).OR..NOT.PRESENT(pcen_in))) & 150 CALL abort_physic(sub, 'Missing lat_in and/or pcen_in (adjust_tropopause=T)', 1) 148 151 149 152 !$OMP MASTER … … 155 158 lPrTrop=NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==NF90_NOERR 156 159 lO3Trop=NF90_INQ_VARID(fID,"tro3_at_tropopause" ,vID)==NF90_NOERR 157 linterp=PRESENT(time_in) 158 IF(linterp) linterp=SIZE(time_in,1)==14 160 linterp=PRESENT(time_in); IF(linterp) linterp=ntim_in==14 159 161 IF(linterp) THEN 160 162 ALLOCATE(v1m(nbp_lon,nbp_lat,nlev_in,n_var)) … … 165 167 END IF 166 168 !--- INITIAL INDEX: LOCATE A LAYER WELL ABOVE TROPOPAUSE (50hPa) 167 IF(PRESENT(pcen_in)) THEN 168 itrp0=locate(pcen_in,prtrop) 169 ELSE 170 itrp0=locate(SQRT(pint_in(1:nlev_in)*pint_in(2:nlev_in+1)),prtrop) 171 END IF 172 IF(lPrSurf) WRITE(*,*)' >> Using GROUND PRESSURE from input O3 forcing file.' 173 IF(linterp) WRITE(*,*)' >> Monthly O3 files => ONLINE TIME INTERPOLATION.' 174 WRITE(*,*)' >> o3 forcing file tropopause location uses:' 175 IF(lPrTrop) THEN; WRITE(*,*)' PRESSURE AT TROPOPAUSE from file.' 176 ELSE IF(lO3Trop) THEN; WRITE(*,*)' O3 CONCENTRATION AT TROPOPAUSE from file.' 177 ELSE; WRITE(*,*)' PARAMETRIZATION of O3 concentration at tropopause.' 169 IF(lAdjTro) itrp0=locate(pcen_in,prtrop) 170 IF(lPrSurf) WRITE(lunout,*)TRIM(sub)//': Using GROUND PRESSURE from input O3 forcing file.' 171 IF(linterp) WRITE(lunout,*)TRIM(sub)//': Monthly O3 files => ONLINE TIME INTERPOLATION.' 172 IF(lAdjTro) THEN 173 WRITE(lunout,*)TRIM(sub)//': o3 forcing file tropopause location uses:' 174 IF(lPrTrop) THEN 175 WRITE(lunout,*)' PRESSURE AT TROPOPAUSE from file.' 176 ELSE IF(lO3Trop) THEN 177 WRITE(lunout,*)' O3 CONCENTRATION AT TROPOPAUSE from file.' 178 ELSE 179 WRITE(lunout,*)' PARAMETRIZED O3 concentration at tropopause.' 180 END IF 178 181 END IF 179 182 END IF 180 183 181 184 !=== UPDATE (ALWAYS FOR DAILY FILES, EACH MONTH FOR MONTHLY FILES) 182 IF(.NOT.linterp.OR.(linterp.AND.(lfirst.OR.julien>=time_in(irec+1)))) &183 CALL update_fields() 185 CALL update_fields() 186 184 187 185 188 !=== TIME INTERPOLATION FOR MONTHLY INPUT FILES 186 189 IF(linterp) THEN 187 WRITE( *,'(3(a,f7.3))')' >> Interpolating O3 at julian day ',julien,&188 & ' from forcingfields at times ',time_in(irec),' and ', time_in(irec+1)190 WRITE(lunout,'(3(a,f7.3))')TRIM(sub)//': Interpolating O3 at julian day '& 191 &,julien,' from fields at times ',time_in(irec),' and ', time_in(irec+1) 189 192 al=(time_in(irec+1)-julien)/(time_in(irec+1)-time_in(irec)) 190 193 v1=al*v1m+(1.-al)*v1p … … 195 198 END IF 196 199 !$OMP END MASTER 200 IF(lfirst) THEN 201 lfirst=.FALSE. ; CALL bcast(lfirst) ; CALL bcast(lPrSurf) 202 CALL bcast(lPrTrop); CALL bcast(lO3Trop); CALL bcast(linterp) 203 IF(lAdjTro) CALL bcast(itrp0) 204 END IF 197 205 CALL scatter2d(v1,v2) 198 CALL bcast(itrp0)199 206 !--- No "ps" in input file => assumed to be equal to current LMDZ ground press 200 207 IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=Pint_ou(:,1); END IF … … 203 210 204 211 !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing 205 IF(.NOT.l trop) THEN212 IF(.NOT.lAdjTro) THEN 206 213 DO i=1,klon 207 214 CALL regr_conserv(1,v2(i,:,:) , Pint_in , Pint_ou(i,nbp_lev+1:1:-1), & … … 245 252 !------------------------------------------------------------------------------- 246 253 IF(.NOT.linterp) THEN !=== DAILY FILES: NO TIME INTERPOLATION 247 WRITE(*,*)' >> Updating Ozone forcing field: read from file.' 248 irec=INT(julien)+1 !--- irec is just the julian day number 254 WRITE(lunout,*)TRIM(sub)//': Updating Ozone forcing field: read from file.' 255 irec=MIN(INT(julien)+1,ntim_in) !--- irec is just the julian day number 256 !--- MIN -> Security in the unlikely case of roundup errors. 249 257 CALL get_3Dfields(v1) !--- Read ozone field(s) 250 IF(l trop) THEN!--- Additional files for fields strain258 IF(lAdjTro) THEN !--- Additional files for fields strain 251 259 IF(lPrSurf) CALL get_2Dfield(ps1,"ps") 252 260 IF(lPrTrop) CALL get_2Dfield(pt1,"tropopause_air_pressure") … … 254 262 END IF 255 263 ELSE !=== MONTHLY FILES: GET 2 NEAREST RECS 256 WRITE(*,*)' >> Refreshing adjacent Ozone forcing fields.' 264 IF(.NOT.lfirst.AND.julien<time_in(irec+1)) RETURN 265 WRITE(lunout,*)TRIM(sub)//': Refreshing adjacent Ozone forcing fields.' 257 266 IF(lfirst) THEN !=== READ EARLIEST TIME FIELDS 258 267 irec=locate(time_in,julien) !--- Need to locate surrounding times 259 268 CALL get_3Dfields(v1m) !--- Read ozone field(s) 260 IF(l trop) THEN!--- Additional files for fields strain269 IF(lAdjTro) THEN !--- Additional files for fields strain 261 270 IF(lPrSurf) CALL get_2Dfield(psm,"ps") 262 271 IF(lPrTrop) CALL get_2Dfield(ptm,"tropopause_air_pressure") … … 265 274 ELSE !=== SHIFT FIELDS 266 275 v1m=v1p !--- Ozone fields 267 IF(l trop) THEN!--- Additional files for fields strain276 IF(lAdjTro) THEN !--- Additional files for fields strain 268 277 IF(lPrSurf) psm=psp !--- Surface pressure 269 278 IF(lPrTrop) ptm=ptp !--- Tropopause pressure … … 273 282 irec=irec+1 274 283 CALL get_3Dfields(v1p) !--- Read ozone field(s) 275 IF(l trop) THEN!--- Additional files for fields strain284 IF(lAdjTro) THEN !--- Additional files for fields strain 276 285 IF(lPrSurf) CALL get_2Dfield(psp,"ps") 277 286 IF(lPrTrop) CALL get_2Dfield(ptp,"tropopause_air_pressure") … … 280 289 IF(lfirst) irec=irec-1 281 290 END IF 282 lfirst=.FALSE.283 CALL bcast(lfirst)284 291 285 292 END SUBROUTINE update_fields
Note: See TracChangeset
for help on using the changeset viewer.
