[11589] | 1 | PROGRAM main |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** PROGRAM main *** |
---|
| 4 | !! |
---|
| 5 | !! ** Purpose : horizontal extrapolation ("drowning") of ECMWF dataset |
---|
| 6 | !! |
---|
| 7 | !! |
---|
| 8 | !!====================================================================== |
---|
| 9 | !! History : 2016-10 (F. Lemarié) Original code largely inspired by SOSIE (L. Brodeau & J.M. Molines) |
---|
| 10 | !! |
---|
| 11 | !!---------------------------------------------------------------------- |
---|
| 12 | USE module_io ! I/O routines |
---|
| 13 | USE module_grid ! compute input and output grids |
---|
| 14 | !! |
---|
| 15 | IMPLICIT NONE |
---|
| 16 | !!---------------------------------------------------------------------- |
---|
| 17 | !! |
---|
| 18 | !! |
---|
| 19 | !! |
---|
| 20 | !!---------------------------------------------------------------------- |
---|
| 21 | ! |
---|
| 22 | INTEGER :: ji,jj,jk,kt,kv,jjp1,jjm1,jip1,jim1 |
---|
| 23 | INTEGER :: jpka,jpvar ! number of vertical levels for input and target grids |
---|
| 24 | INTEGER :: jpi , jpj ! number of grid points in x and y directions |
---|
| 25 | INTEGER :: jptime,ctrl,niter,status,ncid |
---|
| 26 | INTEGER :: ioerr |
---|
| 27 | INTEGER, PARAMETER :: stdout = 6 |
---|
| 28 | INTEGER, PARAMETER :: max_iter = 50 |
---|
| 29 | !! |
---|
| 30 | REAL(8) :: cff,cff1,cff2 |
---|
| 31 | !! |
---|
| 32 | REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: varin |
---|
| 33 | REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: tmask, tmask2, ff_t |
---|
| 34 | REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: mask_coast, maskv ! land-sea mask |
---|
| 35 | REAL(8), ALLOCATABLE, DIMENSION(:,: ) :: varold,varnew |
---|
| 36 | REAL(8), ALLOCATABLE, DIMENSION(: ) :: tmp1d, lat |
---|
| 37 | !! |
---|
| 38 | CHARACTER(len=500) :: file_u,file_v,file_hpg,file_geos ! ERAi files containing wind components |
---|
| 39 | CHARACTER(len=500) :: file_t,file_q, file_m ! ERAi files containing tracers and mask |
---|
| 40 | CHARACTER(len=500) :: file_z,file_p,cn_dir ! ERAi files containing surface geopot and pressure |
---|
| 41 | CHARACTER(len=500) :: abl_file, grd_file, drwn_file, out_file |
---|
| 42 | CHARACTER(len=500) :: namelistf |
---|
| 43 | CHARACTER(len=500) :: argument |
---|
| 44 | CHARACTER(len= 20), DIMENSION(4) :: dimnames |
---|
| 45 | CHARACTER(:), ALLOCATABLE, DIMENSION(:) :: varnames |
---|
| 46 | CHARACTER(6) :: mask_var ! name of mask variable in file_m file |
---|
| 47 | CHARACTER(6) :: var_name |
---|
| 48 | !! |
---|
| 49 | LOGICAL :: ln_read_zsurf ! read surface geopotential or not |
---|
| 50 | LOGICAL :: ln_read_mask ! read land-sea mask or not |
---|
| 51 | LOGICAL :: ln_perio_latbc ! use periodic BC along the domain latitudinal edges (for global data) or use zero-gradient BC (for regional data) |
---|
| 52 | LOGICAL :: ln_c1d ! output only a single column in output file |
---|
| 53 | LOGICAL :: ln_hpg_frc ! compute horizontal pressure gradient |
---|
| 54 | LOGICAL :: ln_geo_wnd ! compute goestrophic wind components |
---|
| 55 | LOGICAL :: ln_slp_smth ! apply gibbs oscillation filetring on mean sea level pressure |
---|
| 56 | LOGICAL :: ln_drw_smth ! apply gibbs oscillation filetring on mean sea level pressure |
---|
| 57 | LOGICAL :: ln_slp_log ! log(sea-level pressure) or sea-level pressure |
---|
| 58 | LOGICAL :: ln_lsm_land ! if T mask is 1 over land and 0 over ocean if F it is the other way around |
---|
| 59 | INTEGER :: ptemp_method ! way to compute potential temperature |
---|
| 60 | !! |
---|
| 61 | REAL(8), PARAMETER :: omega = 7.292116e-05 |
---|
| 62 | REAL(8), PARAMETER :: rad = 3.141592653589793 / 180. |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | !!--------------------------------------------------------------------- |
---|
| 66 | !! List of variables read in the namelist file |
---|
| 67 | NAMELIST/nml_out/ grd_file, abl_file, drwn_file, var_name |
---|
| 68 | NAMELIST/nml_opt/ ptemp_method, ln_slp_log, ln_slp_smth, ln_read_mask, ln_perio_latbc, & |
---|
| 69 | & ln_hpg_frc, ln_geo_wnd, ln_c1d, ln_read_zsurf, ln_lsm_land, ln_drw_smth |
---|
| 70 | NAMELIST/nml_fld/ cn_dir, file_u, file_v, file_t, & |
---|
| 71 | & file_q, file_z, file_p, file_hpg, file_geos, & |
---|
| 72 | & file_m, mask_var |
---|
| 73 | !! |
---|
| 74 | !! get the namelist file name |
---|
| 75 | CALL get_command_argument(1, argument, ctrl, status) |
---|
| 76 | ! |
---|
| 77 | SELECT CASE(status) |
---|
| 78 | CASE(0) |
---|
| 79 | namelistf = trim(argument) |
---|
| 80 | CASE(-1) |
---|
| 81 | WRITE(stdout,*) "### Error: file name too long" |
---|
| 82 | STOP |
---|
| 83 | CASE DEFAULT |
---|
| 84 | namelistf = 'namelist_abl_tools' |
---|
| 85 | END SELECT |
---|
| 86 | !!--------------------------------------------------------------------- |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | !!--------------------------------------------------------------------- |
---|
| 90 | !! read namelist variables |
---|
| 91 | ctrl = 0 |
---|
| 92 | OPEN(50, file=namelistf, status='old', form='formatted', access='sequential', iostat=ioerr) |
---|
| 93 | IF (ioerr /= 0) ctrl = ctrl + 1 |
---|
| 94 | READ(50,nml_opt, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1 |
---|
| 95 | READ(50,nml_fld, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1 |
---|
| 96 | READ(50,nml_out, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1 |
---|
| 97 | |
---|
| 98 | IF (ctrl > 0) then |
---|
| 99 | WRITE(stdout,*) "### E R R O R while reading namelist file '",trim(namelistf),"'" |
---|
| 100 | WRITE(stdout,*) " ctrl = ",ctrl |
---|
| 101 | STOP |
---|
| 102 | ELSE |
---|
| 103 | WRITE(stdout,*) " Namelist file ",trim(namelistf)," OK " |
---|
| 104 | END IF |
---|
| 105 | !!------------------------------------------------------------------------------------- |
---|
| 106 | !! list of variables to treat |
---|
| 107 | !! |
---|
| 108 | !! get the variable name |
---|
| 109 | CALL get_command_argument( 2, argument, ctrl, status) |
---|
| 110 | SELECT CASE(status) |
---|
| 111 | CASE(0) |
---|
| 112 | var_name = trim(argument) |
---|
| 113 | CASE(-1) |
---|
| 114 | WRITE(stdout,*) "### Error: file name too long" |
---|
| 115 | STOP |
---|
| 116 | END SELECT |
---|
| 117 | WRITE(stdout,*) "var_name: ", trim(var_name), Len_Trim(var_name) |
---|
| 118 | ! |
---|
| 119 | IF (Len_Trim(var_name) == 0) THEN |
---|
| 120 | jpvar = 7 |
---|
| 121 | allocate(character(len=20) :: varnames(jpvar) ) |
---|
| 122 | varnames(1:5) = [character(len=4) :: 'slp', 'uwnd', 'vwnd', 'tair', 'humi' ] |
---|
| 123 | IF (ln_hpg_frc) varnames(6:7) = [character(len=4) :: 'uhpg', 'vhpg' ] |
---|
| 124 | IF (ln_geo_wnd) varnames(6:7) = [character(len=4) :: 'ugeo', 'vgeo' ] |
---|
| 125 | IF( .NOT. Var_Existence( varnames(5), trim(cn_dir)//'/'//trim(abl_file) ) ) jpvar = jpvar - 1 |
---|
| 126 | IF( .NOT. Var_Existence( varnames(6), trim(cn_dir)//'/'//trim(abl_file) ) ) jpvar = jpvar - 1 |
---|
| 127 | ELSE |
---|
| 128 | jpvar = 1 |
---|
| 129 | allocate(character(len=20) :: varnames(jpvar) ) |
---|
| 130 | varnames = var_name |
---|
| 131 | abl_file = trim(var_name)//'_'//trim( abl_file) |
---|
| 132 | drwn_file = trim(var_name)//'_'//trim(drwn_file) |
---|
| 133 | END IF |
---|
| 134 | WRITE(stdout,*) "Number of variables to treat : ",jpvar |
---|
| 135 | !!------------------------------------------------------------------------------------- |
---|
| 136 | !! read the dimensions for the input files |
---|
| 137 | CALL Read_Ncdf_dim ( 'jpka' , trim(cn_dir)//'/'//trim(abl_file), jpka ) |
---|
| 138 | CALL Read_Ncdf_dim ( 'time' , trim(cn_dir)//'/'//trim(abl_file), jptime ) |
---|
| 139 | CALL Read_Ncdf_dim ( 'lon' , trim(cn_dir)//'/'//trim(abl_file), jpi ) |
---|
| 140 | CALL Read_Ncdf_dim ( 'lat' , trim(cn_dir)//'/'//trim(abl_file), jpj ) |
---|
| 141 | WRITE(stdout,*) "jpka, jptime, jpi, jpj: ", jpka, jptime, jpi, jpj |
---|
| 142 | ! |
---|
| 143 | !!--------------------------------------------------------------------- |
---|
| 144 | |
---|
| 145 | |
---|
| 146 | !!--------------------------------------------------------------------- |
---|
| 147 | !! allocate arrays |
---|
| 148 | ALLOCATE(varin (1:jpi,1:jpj,1:jpka,1)) |
---|
| 149 | ALLOCATE(tmask (1:jpi,1:jpj )) |
---|
| 150 | ALLOCATE(tmask2(1:jpi,1:jpj )) |
---|
| 151 | ALLOCATE(maskv (1:jpi,1:jpj )) |
---|
| 152 | ALLOCATE(mask_coast (1:jpi,1:jpj )) |
---|
| 153 | ALLOCATE(varold (1:jpi,1:jpj )) |
---|
| 154 | ALLOCATE(varnew (1:jpi,1:jpj )) |
---|
| 155 | |
---|
| 156 | !!--------------------------------------------------------------------- |
---|
| 157 | !! Read the mask and remove some closed seas |
---|
| 158 | IF (ln_read_mask) THEN |
---|
| 159 | CALL init_atm_mask(jpi,jpj,trim(cn_dir)//'/'//trim(file_m),trim(mask_var),ln_lsm_land, tmask ) |
---|
| 160 | CALL Read_Ncdf_var ( 'tmask' , trim(cn_dir)//'/'//trim(abl_file), tmask2(:,:) ) |
---|
| 161 | ELSE |
---|
| 162 | tmask(:,:) = 1. |
---|
| 163 | tmask2(:,:) = 1. |
---|
| 164 | END IF |
---|
| 165 | !! |
---|
| 166 | |
---|
| 167 | |
---|
| 168 | |
---|
| 169 | !!--------------------------------------------------------------------- |
---|
| 170 | !! create output file |
---|
| 171 | !! |
---|
| 172 | out_file = trim(cn_dir)//'/'//trim(drwn_file) |
---|
| 173 | CALL Init_output_File ( jpi, jpj, jpka-1, trim(cn_dir)//'/'//trim(abl_file), out_file, tmask(:,:) ) |
---|
| 174 | |
---|
| 175 | !!--------------------------------------------------------------------- |
---|
| 176 | !! Initialize the name of the dimensions for the result of the drowning |
---|
| 177 | !! |
---|
| 178 | dimnames(1) = 'lon' |
---|
| 179 | dimnames(2) = 'lat' |
---|
| 180 | dimnames(3) = 'jpka' |
---|
| 181 | dimnames(4) = 'time' |
---|
| 182 | |
---|
| 183 | CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask2, 'float' ) |
---|
| 184 | |
---|
| 185 | !!--------------------------------------------------------------------- |
---|
| 186 | ! Read time variable |
---|
| 187 | ALLOCATE(tmp1d (1:jptime)) |
---|
| 188 | ALLOCATE(lat (1:jpj)) |
---|
| 189 | CALL Read_Ncdf_var ( 'time', trim(cn_dir)//'/'//trim(file_t), tmp1d ) |
---|
| 190 | CALL Read_Ncdf_var ( 'lat' , trim(cn_dir)//'/'//trim(file_t), lat ) !<-- latitude |
---|
| 191 | !!--------------------------------------------------------------------- |
---|
| 192 | |
---|
| 193 | ! force drowning in the equatorial band with geo wind |
---|
| 194 | IF (ln_geo_wnd) THEN |
---|
| 195 | ALLOCATE( ff_t(1:jpi,1:jpj) ) |
---|
| 196 | DO jj = 1, jpj |
---|
| 197 | DO ji = 1, jpi |
---|
| 198 | ff_t(ji,jj) = 2. * omega * SIN( rad * lat(jj) ) |
---|
| 199 | IF (abs(ff_t(ji,jj)) < 2.5e-5) tmask2(ji,jj) = 0. |
---|
| 200 | END DO |
---|
| 201 | END DO |
---|
| 202 | END IF |
---|
| 203 | |
---|
| 204 | !=========== |
---|
| 205 | DO kv = 1,jpvar |
---|
| 206 | !=========== |
---|
| 207 | DO kt = 1,jptime |
---|
| 208 | !=========== |
---|
| 209 | |
---|
| 210 | IF( kv == 1 ) THEN |
---|
| 211 | CALL Write_Ncdf_var( 'time', dimnames(4:4), trim(out_file), tmp1d(kt:kt), kt, 'double' ) |
---|
| 212 | ENDIF |
---|
| 213 | |
---|
| 214 | ! Read variable to treat |
---|
| 215 | print*,'Treat variable ',trim(varnames(kv)),' at time level ',kt |
---|
| 216 | IF (trim(varnames(kv)).EQ."slp") THEN |
---|
| 217 | CALL Read_Ncdf_var ( trim(varnames(kv)), trim(cn_dir)//'/'//trim(file_hpg), varin(:,:,1,1), kt ) |
---|
| 218 | ELSE |
---|
| 219 | CALL Read_Ncdf_var ( trim(varnames(kv)), trim(cn_dir)//'/'//trim(abl_file), varin, kt ) |
---|
| 220 | END IF |
---|
| 221 | |
---|
| 222 | |
---|
| 223 | !=========== |
---|
| 224 | DO jk = 1,jpka |
---|
| 225 | !=========== |
---|
| 226 | |
---|
| 227 | IF ( (trim(varnames(kv))=="uwnd").OR.(trim(varnames(kv))=="vwnd").OR. & |
---|
| 228 | & (trim(varnames(kv))=="ugeo").OR.(trim(varnames(kv))=="vgeo").OR. & |
---|
| 229 | & (trim(varnames(kv))=="uhpg").OR.(trim(varnames(kv))=="vhpg") ) THEN |
---|
| 230 | maskv (1:jpi,1:jpj ) = tmask2(1:jpi,1:jpj ) |
---|
| 231 | ELSE |
---|
| 232 | maskv (1:jpi,1:jpj ) = tmask(1:jpi,1:jpj ) |
---|
| 233 | END IF |
---|
| 234 | |
---|
| 235 | |
---|
| 236 | varold (1:jpi,1:jpj) = varin(1:jpi,1:jpj,jk,1) |
---|
| 237 | varnew (1:jpi,1:jpj) = varin(1:jpi,1:jpj,jk,1) |
---|
| 238 | mask_coast(1:jpi,1:jpj) = 0. |
---|
| 239 | |
---|
| 240 | !$$$$$$$$$$$$$$$$ |
---|
| 241 | DO niter = 1, max_iter |
---|
| 242 | !$$$$$$$$$$$$$$$$ |
---|
| 243 | |
---|
| 244 | varold(1:jpi,1:jpj) = varnew (1:jpi,1:jpj) |
---|
| 245 | |
---|
| 246 | IF ( .NOT. (ANY(maskv == 0)) ) THEN |
---|
| 247 | EXIT |
---|
| 248 | END IF |
---|
| 249 | |
---|
| 250 | !+++++++++++++++++++++++++++++++++++ |
---|
| 251 | ! Build mask_coast |
---|
| 252 | DO jj = 1,jpj |
---|
| 253 | DO ji = 1,jpi |
---|
| 254 | |
---|
| 255 | IF (ln_perio_latbc) THEN |
---|
| 256 | jip1 = ji + 1; if(ji==jpi) jip1 = 1 |
---|
| 257 | jim1 = ji - 1; if(ji== 1) jim1 = jpi |
---|
| 258 | ELSE |
---|
| 259 | jip1 = ji + 1; if(ji==jpi) jip1 = jpi-1 |
---|
| 260 | jim1 = ji - 1; if(ji== 1) jim1 = 2 |
---|
| 261 | END IF |
---|
| 262 | jjp1 = jj + 1; if(jj==jpj) jjp1 = jpj-1 |
---|
| 263 | jjm1 = jj - 1; if(jj== 1) jjm1 = 2 |
---|
| 264 | |
---|
| 265 | cff = (1.-maskv(ji,jj)) |
---|
| 266 | cff1 = maskv(ji,jjp1)+maskv(jip1,jj)+maskv(jim1,jj)+maskv(ji,jjm1) |
---|
| 267 | mask_coast(ji,jj) = cff * cff1 |
---|
| 268 | IF( mask_coast(ji,jj) > 0. ) mask_coast(ji,jj) = 1. |
---|
| 269 | |
---|
| 270 | END DO |
---|
| 271 | END DO |
---|
| 272 | !+++++++++++++++++++++++++++++++++++ |
---|
| 273 | |
---|
| 274 | DO jj=1,jpj |
---|
| 275 | DO ji=1,jpi |
---|
| 276 | |
---|
| 277 | IF ( mask_coast(ji,jj) == 1. ) THEN |
---|
| 278 | |
---|
| 279 | IF (ln_perio_latbc) THEN |
---|
| 280 | jip1 = ji + 1; if(ji==jpi) jip1 = 1 |
---|
| 281 | jim1 = ji - 1; if(ji== 1) jim1 = jpi |
---|
| 282 | ELSE |
---|
| 283 | jip1 = ji + 1; if(ji==jpi) jip1 = jpi-1 |
---|
| 284 | jim1 = ji - 1; if(ji== 1) jim1 = 2 |
---|
| 285 | END IF |
---|
| 286 | jjp1 = jj + 1; if(jj==jpj) jjp1 = jpj-1 |
---|
| 287 | jjm1 = jj - 1; if(jj== 1) jjm1 = 2 |
---|
| 288 | |
---|
| 289 | cff = maskv(jim1,jjm1)+maskv(jim1,jj)+maskv(jim1,jjp1) & |
---|
| 290 | & + maskv(ji ,jjm1)+ maskv(ji ,jjp1) & |
---|
| 291 | & + maskv(jip1,jjm1)+maskv(jip1,jj)+maskv(jip1,jjp1) |
---|
| 292 | |
---|
| 293 | varnew(ji,jj) = (1./cff)*( & |
---|
| 294 | & varold(jim1,jjm1)*maskv(jim1,jjm1) & |
---|
| 295 | & + varold(jim1,jj )*maskv(jim1,jj ) & |
---|
| 296 | & + varold(jim1,jjp1)*maskv(jim1,jjp1) & |
---|
| 297 | & + varold(ji ,jjm1)*maskv(ji ,jjm1) & |
---|
| 298 | & + varold(ji ,jjp1)*maskv(ji ,jjp1) & |
---|
| 299 | & + varold(jip1,jjm1)*maskv(jip1,jjm1) & |
---|
| 300 | & + varold(jip1,jj )*maskv(jip1,jj ) & |
---|
| 301 | & + varold(jip1,jjp1)*maskv(jip1,jjp1) ) |
---|
| 302 | |
---|
| 303 | END IF |
---|
| 304 | |
---|
| 305 | END DO |
---|
| 306 | END DO |
---|
| 307 | !+++++++++++++++++++++++++++++++++++ |
---|
| 308 | maskv(1:jpi,1:jpj) = maskv(1:jpi,1:jpj) + mask_coast(1:jpi,1:jpj) |
---|
| 309 | !+++++++++++++++++++++++++++++++++++ |
---|
| 310 | !$$$$$$$$$$$$$$$$ |
---|
| 311 | END DO |
---|
| 312 | !$$$$$$$$$$$$$$$$ |
---|
| 313 | |
---|
| 314 | ! SMOOTHING AFTER DROWNING |
---|
| 315 | IF (ln_drw_smth) THEN |
---|
| 316 | maskv (1:jpi,1:jpj ) = 1. |
---|
| 317 | ELSE |
---|
| 318 | !!only over continent |
---|
| 319 | maskv (1:jpi,1:jpj ) = 1. - tmask(1:jpi,1:jpj ) |
---|
| 320 | END IF |
---|
| 321 | |
---|
| 322 | CALL smooth_field( jpi, jpj, varnew, maskv, 3 ) |
---|
| 323 | varin(1:jpi,1:jpj,jk,1) = varnew(1:jpi,1:jpj) |
---|
| 324 | |
---|
| 325 | IF (trim(varnames(kv)).EQ."slp") EXIT |
---|
| 326 | |
---|
| 327 | !=========== |
---|
| 328 | END DO ! jk |
---|
| 329 | !=========== |
---|
| 330 | |
---|
| 331 | IF (trim(varnames(kv)).EQ."slp") THEN |
---|
| 332 | CALL Write_Ncdf_var ( trim(varnames(kv)), dimnames(1:4), trim(cn_dir)//'/'//drwn_file, varin(:,:,1,1), kt, 'float' ) |
---|
| 333 | ELSE |
---|
| 334 | CALL Write_Ncdf_var ( trim(varnames(kv)), dimnames(1:4), trim(cn_dir)//'/'//drwn_file, varin(:,:,:,:), kt, 'float' ) |
---|
| 335 | END IF |
---|
| 336 | |
---|
| 337 | !=========== |
---|
| 338 | END DO ! time |
---|
| 339 | !=========== |
---|
| 340 | |
---|
| 341 | !=========== |
---|
| 342 | END DO ! variable |
---|
| 343 | !=========== |
---|
| 344 | |
---|
| 345 | STOP |
---|
| 346 | ! |
---|
| 347 | END PROGRAM main |
---|