[8] | 1 | PROGRAM driver |
---|
| 2 | !- $Header: /home/ssipsl/CVSREP/ORCHIDEE_OL/dim2_driver.f90,v 1.26 2010/04/22 15:48:59 ssipsl Exp $ |
---|
| 3 | !- IPSL (2006) |
---|
| 4 | !- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
| 5 | !--------------------------------------------------------------------- |
---|
| 6 | !- This PROGRAM is the driver of the dim2 version of SECHIBA. It's |
---|
| 7 | !- main use is for PILPS type experiments. To run it you need to have |
---|
| 8 | !- the following software : |
---|
| 9 | !- - SECHIBA |
---|
| 10 | !- - ioipsl |
---|
| 11 | !- - netcdf |
---|
| 12 | !- - F90 compiler |
---|
| 13 | !- - tk (optional but very useful for configuring the model) |
---|
| 14 | !- Obviously also the associated Makefiles. |
---|
| 15 | !- |
---|
| 16 | !- The forcing data needs to be in netCDF format and should |
---|
| 17 | !- contain the following variables : |
---|
| 18 | !- - Incoming SW radiation |
---|
| 19 | !- - Incoming LW radiation |
---|
| 20 | !- - Precipitation |
---|
| 21 | !- - Air temperature at a reference level |
---|
| 22 | !- - Air humidity at the same level |
---|
| 23 | !- - wind at the same level |
---|
| 24 | !- - surface pressure |
---|
| 25 | !- |
---|
| 26 | !- Once you have all this and compiled the model you can configure it. |
---|
| 27 | !- Type make config and set all the values presented in the menu. |
---|
| 28 | !- This tool will create a run.def which will be used by the model. |
---|
| 29 | !- The run.def can also be edited by hand but it is more tedious. |
---|
| 30 | !- Once this run.def is created the model is ready to run. |
---|
| 31 | !--------------------------------------------------------------------- |
---|
| 32 | USE netcdf |
---|
| 33 | USE ioipsl |
---|
| 34 | USE grid |
---|
| 35 | USE intersurf, ONLY : intersurf_main |
---|
| 36 | USE constantes |
---|
| 37 | USE constantes_veg |
---|
| 38 | USE readdim2 |
---|
| 39 | USE parallel |
---|
| 40 | USE timer |
---|
| 41 | !- |
---|
| 42 | IMPLICIT NONE |
---|
| 43 | !- |
---|
| 44 | INTEGER :: iim, jjm, llm |
---|
| 45 | INTEGER :: im, jm, lm, tm, is, force_id, itest, jtest |
---|
| 46 | REAL :: dt, dt_force, dt_rest, date0, date0_rest |
---|
| 47 | REAL :: zlflu |
---|
| 48 | REAL :: alpha |
---|
| 49 | !- |
---|
| 50 | REAL, ALLOCATABLE, DIMENSION(:,:) :: & |
---|
| 51 | & swdown, precip_rain, precip_snow, tair_obs, u, v, & |
---|
| 52 | & qair_obs, pb, lwdown, & |
---|
| 53 | & eair_obs, zlev_vec, zlevuv_vec, relax |
---|
| 54 | !- Variables which are forcings for SECHIBA |
---|
| 55 | REAL, ALLOCATABLE, DIMENSION(:,:) :: & |
---|
| 56 | & petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, & |
---|
| 57 | & for_u, for_v, for_swnet, for_swdown, for_lwdown, & |
---|
| 58 | & for_psurf, for_qair, for_tair, for_eair, & |
---|
| 59 | & for_ccanopy, for_rau |
---|
| 60 | !!$, tmp_eair, tmp_tair, & |
---|
| 61 | !!$ & tmp_qair, tmp_pb |
---|
| 62 | !- |
---|
| 63 | REAL, ALLOCATABLE, DIMENSION(:,:) :: & |
---|
| 64 | & for_contfrac, old_zlev, old_qair, old_eair, tsol_rad, vevapp, & |
---|
| 65 | & temp_sol_NEW, temp_sol_old, qsurf, dtdt, coastalflow, riverflow, & |
---|
| 66 | & fluxsens, fluxlat, emis, z0 |
---|
| 67 | !!$, epot_sol |
---|
| 68 | !- |
---|
| 69 | INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: for_neighbours |
---|
| 70 | !- |
---|
| 71 | REAL, ALLOCATABLE, DIMENSION(:,:,:) :: for_resolution |
---|
| 72 | !- |
---|
| 73 | REAL, ALLOCATABLE, DIMENSION(:,:,:) :: albedo |
---|
| 74 | REAL, ALLOCATABLE, DIMENSION(:,:) :: albedo_vis |
---|
| 75 | REAL, ALLOCATABLE, DIMENSION(:,:) :: albedo_nir |
---|
| 76 | !- |
---|
| 77 | INTEGER, ALLOCATABLE, DIMENSION(:) :: kindex |
---|
| 78 | REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat |
---|
| 79 | !- |
---|
| 80 | REAL :: old_tair |
---|
| 81 | REAL :: atmco2 |
---|
| 82 | INTEGER :: nbindex |
---|
| 83 | REAL,DIMENSION(:), ALLOCATABLE :: lev, levuv |
---|
| 84 | REAL :: julian, ss |
---|
| 85 | INTEGER :: yy, mm, dd |
---|
| 86 | !- |
---|
| 87 | LOGICAL :: relaxation, lower_wind |
---|
| 88 | !- |
---|
| 89 | CHARACTER(LEN=80) :: filename, restname_in, restname_out |
---|
| 90 | CHARACTER(LEN=30) :: time_str, var_name |
---|
| 91 | !- |
---|
| 92 | INTEGER :: it, istp, istp_old, rest_id, it_force |
---|
| 93 | !- |
---|
| 94 | INTEGER :: split, split_start, nb_spread, for_offset |
---|
| 95 | INTEGER :: itau_dep, itau_dep_rest, itau_fin, itau_skip, itau_len |
---|
| 96 | |
---|
| 97 | INTEGER,DIMENSION(2) :: ml |
---|
| 98 | !- |
---|
| 99 | LOGICAL :: last_CALL, first_CALL, debug |
---|
| 100 | LOGICAL :: no_inter, inter_lin, netrad_cons |
---|
| 101 | !- |
---|
| 102 | ! to check variables passed to intersurf |
---|
| 103 | INTEGER :: ik |
---|
| 104 | INTEGER :: i,j |
---|
| 105 | LOGICAL, PARAMETER :: longprint = .FALSE. |
---|
| 106 | INTEGER :: ierr |
---|
| 107 | |
---|
| 108 | REAL, ALLOCATABLE, DIMENSION(:,:) :: & |
---|
| 109 | & fluxsens_g,vevapp_g,old_zlev_g,old_qair_g,old_eair_g,for_rau_g, & |
---|
| 110 | & petAcoef_g, petBcoef_g,peqAcoef_g,peqBcoef_g,albedo_g,z0_g |
---|
| 111 | LOGICAL :: Flag |
---|
| 112 | |
---|
| 113 | REAL :: fill_init |
---|
| 114 | |
---|
| 115 | fill_init=REAL(nf90_fill_real,r_std) |
---|
| 116 | CALL ioconf_expval(val_exp) |
---|
| 117 | !- |
---|
| 118 | ! Init parallelism |
---|
| 119 | |
---|
| 120 | CALL init_para(.FALSE.) |
---|
| 121 | |
---|
| 122 | ! driver only for process root |
---|
| 123 | |
---|
| 124 | !--------------------------------------------------------------------- |
---|
| 125 | !- |
---|
| 126 | ! set debug to have more information |
---|
| 127 | !- |
---|
| 128 | !Config Key = DEBUG_INFO |
---|
| 129 | !Config Desc = Flag for debug information |
---|
| 130 | !Config Def = n |
---|
| 131 | !Config Help = This option allows to switch on the output of debug |
---|
| 132 | !Config information without recompiling the code. |
---|
| 133 | !- |
---|
| 134 | debug = .FALSE. |
---|
| 135 | CALL getin_p('DEBUG_INFO',debug) |
---|
| 136 | !===================================================================== |
---|
| 137 | !- 1.0 This section defines the general set-up of the experiment : |
---|
| 138 | !- - Forcing data to be used with its time step |
---|
| 139 | !- - Restart file to be used |
---|
| 140 | !- - The time step that will be used |
---|
| 141 | !- - The starting date of the simulation |
---|
| 142 | !- - Length of integration |
---|
| 143 | !- - Basic flags for SSIPSL |
---|
| 144 | !===================================================================== |
---|
| 145 | !- 1.1 Initialize the driving variables. It essentialy gets the mode |
---|
| 146 | !- used and the size of the driving variables. |
---|
| 147 | !===================================================================== |
---|
| 148 | IF (debug) WRITE(numout,*) 'Reading name of the forcing file' |
---|
| 149 | !- |
---|
| 150 | !Config Key = FORCING_FILE |
---|
| 151 | !Config Desc = Name of file containing the forcing data |
---|
| 152 | !Config Def = islscp_for.nc |
---|
| 153 | !Config Help = This is the name of the file which should be opened |
---|
| 154 | !Config for reading the forcing data of the dim0 model. |
---|
| 155 | !Config The format of the file has to be netCDF and COADS |
---|
| 156 | !Config compliant. |
---|
| 157 | !- |
---|
| 158 | filename='islscp_for.nc' |
---|
| 159 | CALL getin_p('FORCING_FILE',filename) |
---|
| 160 | !- |
---|
| 161 | IF (debug) WRITE(numout,*) 'Opening forcing file' |
---|
| 162 | !- |
---|
| 163 | ! We call flininfo to obtain the dimensions |
---|
| 164 | ! of iim, jjm and others. |
---|
| 165 | ! This information will allow us to allocate all the space needed. |
---|
| 166 | !- |
---|
| 167 | CALL forcing_info & |
---|
| 168 | & (filename, iim, jjm, llm, tm, date0, dt_force, force_id) |
---|
| 169 | !- |
---|
| 170 | WRITE(numout,*) 'Out of flininfo : date0 ', date0, & |
---|
| 171 | 'iim, jjm, llm, tm',iim,jjm,llm,tm,' dt_force ',dt_force |
---|
| 172 | !- |
---|
| 173 | IF (debug) THEN |
---|
| 174 | WRITE(numout,*) 'Allocate memory for the driver :', iim, jjm, llm |
---|
| 175 | ENDIF |
---|
| 176 | !- |
---|
| 177 | ALLOCATE(lev(llm),levuv(llm)) |
---|
| 178 | ALLOCATE & |
---|
| 179 | & (swdown(iim,jjm), precip_rain(iim,jjm), precip_snow(iim,jjm), & |
---|
| 180 | & tair_obs(iim,jjm), u(iim,jjm), v(iim,jjm), qair_obs(iim,jjm), & |
---|
| 181 | & pb(iim,jjm), lwdown(iim,jjm), & |
---|
| 182 | & eair_obs(iim,jjm), zlev_vec(iim,jjm), zlevuv_vec(iim,jjm), relax(iim,jjm)) |
---|
| 183 | !- |
---|
| 184 | ALLOCATE & |
---|
| 185 | & (petAcoef(iim,jjm), peqAcoef(iim,jjm), & |
---|
| 186 | & petBcoef(iim,jjm), peqBcoef(iim,jjm), & |
---|
| 187 | & cdrag(iim,jjm), for_u(iim,jjm), for_v(iim,jjm), & |
---|
| 188 | & for_swnet(iim,jjm), for_swdown(iim,jjm), for_lwdown(iim,jjm), & |
---|
| 189 | & for_psurf(iim,jjm), for_qair(iim,jjm), for_tair(iim,jjm), & |
---|
| 190 | & for_eair(iim,jjm), for_ccanopy(iim,jjm), for_rau(iim,jjm)) |
---|
| 191 | !!$, & |
---|
| 192 | !!$ & tmp_eair(iim,jjm), tmp_tair(iim,jjm), & |
---|
| 193 | !!$ & tmp_qair(iim,jjm), tmp_pb(iim,jjm)) |
---|
| 194 | !- |
---|
| 195 | ALLOCATE & |
---|
| 196 | & (for_contfrac(iim,jjm), for_neighbours(iim,jjm,8), for_resolution(iim,jjm,2), & |
---|
| 197 | & old_zlev(iim,jjm), old_qair(iim,jjm), old_eair(iim,jjm), & |
---|
| 198 | & tsol_rad(iim,jjm), vevapp(iim,jjm), & |
---|
| 199 | & temp_sol_NEW(iim,jjm), temp_sol_old(iim,jjm), & |
---|
| 200 | & dtdt(iim,jjm), coastalflow(iim,jjm), riverflow(iim,jjm), & |
---|
| 201 | & fluxsens(iim,jjm), fluxlat(iim,jjm), emis(iim,jjm), & |
---|
| 202 | & z0(iim,jjm), qsurf(iim,jjm)) |
---|
| 203 | !!$, epot_sol(iim,jjm) |
---|
| 204 | !- |
---|
| 205 | ALLOCATE(albedo(iim,jjm,2)) |
---|
| 206 | ALLOCATE(albedo_vis(iim,jjm),albedo_nir(iim,jjm)) |
---|
| 207 | !- |
---|
| 208 | ALLOCATE(kindex(iim*jjm)) |
---|
| 209 | ALLOCATE(lon(iim,jjm), lat(iim,jjm)) |
---|
| 210 | !-- |
---|
| 211 | lev(:) = 0.0 |
---|
| 212 | levuv(:) = 0.0 |
---|
| 213 | swdown(:,:) = fill_init |
---|
| 214 | precip_rain(:,:) = 0.0 |
---|
| 215 | precip_snow(:,:) = 0.0 |
---|
| 216 | tair_obs(:,:) = 0.0 |
---|
| 217 | u(:,:) = fill_init |
---|
| 218 | v(:,:) = fill_init |
---|
| 219 | qair_obs(:,:) = fill_init |
---|
| 220 | pb(:,:) = fill_init |
---|
| 221 | lwdown(:,:) = fill_init |
---|
| 222 | eair_obs(:,:) = fill_init |
---|
| 223 | zlev_vec(:,:) = 0.0 |
---|
| 224 | zlevuv_vec(:,:) = 0.0 |
---|
| 225 | relax(:,:) = 0.0 |
---|
| 226 | petAcoef(:,:) = 0.0 |
---|
| 227 | peqAcoef(:,:) = 0.0 |
---|
| 228 | petBcoef(:,:) = 0.0 |
---|
| 229 | peqBcoef(:,:) = 0.0 |
---|
| 230 | cdrag(:,:) = 0.0 |
---|
| 231 | for_u(:,:) = fill_init |
---|
| 232 | for_v(:,:) = fill_init |
---|
| 233 | for_swnet(:,:) = fill_init |
---|
| 234 | for_swdown(:,:) = fill_init |
---|
| 235 | for_lwdown(:,:) = fill_init |
---|
| 236 | for_psurf(:,:) = fill_init |
---|
| 237 | for_qair(:,:) = fill_init |
---|
| 238 | for_tair(:,:) = fill_init |
---|
| 239 | for_eair(:,:) = fill_init |
---|
| 240 | for_ccanopy(:,:) = 0.0 |
---|
| 241 | for_rau(:,:) = fill_init |
---|
| 242 | for_contfrac(:,:) = 0.0 |
---|
| 243 | for_neighbours(:,:,:) = 0 |
---|
| 244 | for_resolution(:,:,:) = 0.0 |
---|
| 245 | old_zlev(:,:) = 0.0 |
---|
| 246 | old_qair(:,:) = 0.0 |
---|
| 247 | old_eair(:,:) = 0.0 |
---|
| 248 | tsol_rad(:,:) = 0.0 |
---|
| 249 | vevapp(:,:) = 0.0 |
---|
| 250 | temp_sol_NEW(:,:) = fill_init |
---|
| 251 | temp_sol_old(:,:) = fill_init |
---|
| 252 | dtdt(:,:) = 0.0 |
---|
| 253 | coastalflow(:,:) = 0.0 |
---|
| 254 | riverflow(:,:) = 0.0 |
---|
| 255 | fluxsens(:,:) = 0.0 |
---|
| 256 | fluxlat(:,:) = 0.0 |
---|
| 257 | emis(:,:) = 0.0 |
---|
| 258 | z0(:,:) = fill_init |
---|
| 259 | qsurf(:,:) = 0.0 |
---|
| 260 | albedo(:,:,:) = fill_init |
---|
| 261 | albedo_vis(:,:) = fill_init |
---|
| 262 | albedo_nir(:,:) = fill_init |
---|
| 263 | kindex(:) = 0 |
---|
| 264 | lon(:,:) = 0.0 |
---|
| 265 | lat(:,:) = 0.0 |
---|
| 266 | !- |
---|
| 267 | ! We need to know the grid. |
---|
| 268 | ! Then we can initialize the restart files, and then we |
---|
| 269 | ! can give read the restart files in the forcing subroutine. |
---|
| 270 | !- |
---|
| 271 | CALL forcing_grid (iim,jjm,llm,lon,lat,lev,levuv,init_f=.FALSE.) |
---|
| 272 | !===================================================================== |
---|
| 273 | !- 1.2 Time step to be used. |
---|
| 274 | !- This is relative to the time step of the forcing data |
---|
| 275 | !===================================================================== |
---|
| 276 | IF ( .NOT. weathergen ) THEN |
---|
| 277 | !Config Key = SPLIT_DT |
---|
| 278 | !Config Desc = splits the timestep imposed by the forcing |
---|
| 279 | !Config If NOT WEATHERGEN |
---|
| 280 | !Config Def = 12 |
---|
| 281 | !Config Help = With this value the time step of the forcing |
---|
| 282 | !Config will be devided. In principle this can be run |
---|
| 283 | !Config in explicit mode but it is strongly suggested |
---|
| 284 | !Config to use the implicit method so that the |
---|
| 285 | !Config atmospheric forcing has a smooth evolution. |
---|
| 286 | !- |
---|
| 287 | split = 12 |
---|
| 288 | CALL getin_p('SPLIT_DT', split) |
---|
| 289 | ELSE |
---|
| 290 | split = 1 |
---|
| 291 | ENDIF |
---|
| 292 | !- |
---|
| 293 | ! The time step which is going to be used is computed |
---|
| 294 | ! from the forcing time step and the user input. |
---|
| 295 | !- |
---|
| 296 | IF (split > 0) THEN |
---|
| 297 | dt = dt_force/split |
---|
| 298 | ELSE |
---|
| 299 | split = 1 |
---|
| 300 | dt = dt_force |
---|
| 301 | ENDIF |
---|
| 302 | !===================================================================== |
---|
| 303 | !- 1.3 Initialize the restart file for the driver |
---|
| 304 | !===================================================================== |
---|
| 305 | !Config Key = RESTART_FILEIN |
---|
| 306 | !Config Desc = Name of restart to READ for initial conditions |
---|
| 307 | !Config Def = NONE |
---|
| 308 | !Config Help = This is the name of the file which will be opened |
---|
| 309 | !Config to extract the initial values of all prognostic |
---|
| 310 | !Config values of the model. This has to be a netCDF file. |
---|
| 311 | !Config Not truly COADS compliant. NONE will mean that |
---|
| 312 | !Config no restart file is to be expected. |
---|
| 313 | !- |
---|
| 314 | restname_in = 'NONE' |
---|
| 315 | CALL getin_p('RESTART_FILEIN',restname_in) |
---|
| 316 | if (debug) WRITE(numout,*) 'INPUT RESTART_FILE : ',TRIM(restname_in) |
---|
| 317 | !- |
---|
| 318 | !Config Key = RESTART_FILEOUT |
---|
| 319 | !Config Desc = Name of restart files to be created by the driver |
---|
| 320 | !Config Def = driver_rest_out.nc |
---|
| 321 | !Config Help = This variable give the name for |
---|
| 322 | !Config the restart files. The restart software within |
---|
| 323 | !Config IOIPSL will add .nc if needed |
---|
| 324 | !- |
---|
| 325 | restname_out = 'driver_rest_out.nc' |
---|
| 326 | CALL getin_p('RESTART_FILEOUT', restname_out) |
---|
| 327 | if (debug) WRITE(numout,*) 'OUTPUT RESTART_FILE : ',TRIM(restname_out) |
---|
| 328 | !- |
---|
| 329 | ! We choose some default values for the chronology |
---|
| 330 | ! of the restart file. |
---|
| 331 | !- |
---|
| 332 | dt_rest = dt |
---|
| 333 | date0_rest = date0 |
---|
| 334 | !- |
---|
| 335 | ! Set default values for the start and end of the simulation |
---|
| 336 | ! in the forcing chronology. |
---|
| 337 | !- |
---|
| 338 | itau_dep = 0 |
---|
| 339 | itau_dep_rest = 0 |
---|
| 340 | itau_fin = tm-1 |
---|
| 341 | !- |
---|
| 342 | CALL gather2D(lon,lon_g) |
---|
| 343 | CALL gather2D(lat,lat_g) |
---|
| 344 | !- |
---|
| 345 | if (debug) WRITE(numout,*) & |
---|
| 346 | & 'Before Restart file initialization : itau_dep, itau_fin, date0, dt', & |
---|
| 347 | & itau_dep, itau_fin, date0, dt |
---|
| 348 | |
---|
| 349 | IF (is_root_prc) THEN |
---|
| 350 | CALL restini & |
---|
| 351 | & (restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, & |
---|
| 352 | & restname_out, itau_dep_rest, date0_rest, dt_rest, rest_id) |
---|
| 353 | !- |
---|
| 354 | IF (debug) WRITE(numout,*) & |
---|
| 355 | & 'For Restart file initialization : itau_dep_rest, date0_rest, dt_rest', & |
---|
| 356 | & itau_dep_rest, date0_rest, dt_rest |
---|
| 357 | !- |
---|
| 358 | ! MM When itau_dep_rest is not equal to itau_dep |
---|
| 359 | IF (itau_dep /= itau_dep_rest) THEN |
---|
| 360 | itau_dep = itau_dep_rest |
---|
| 361 | itau_fin = itau_dep+tm-1 |
---|
| 362 | ENDIF |
---|
| 363 | ENDIF |
---|
| 364 | CALL bcast (itau_dep_rest) |
---|
| 365 | CALL bcast (itau_dep) |
---|
| 366 | CALL bcast (itau_fin) |
---|
| 367 | CALL bcast (date0_rest) |
---|
| 368 | CALL bcast (dt_rest) |
---|
| 369 | !- |
---|
| 370 | WRITE(numout,*) & |
---|
| 371 | & 'Restart file initialized : itau_dep, itau_fin, date0_rest, dt_rest', & |
---|
| 372 | & itau_dep, itau_fin, date0_rest, dt_rest |
---|
| 373 | !===================================================================== |
---|
| 374 | !- 1.4 Here we do the first real reading of the driving. It only |
---|
| 375 | !- gets a few variables. |
---|
| 376 | !===================================================================== |
---|
| 377 | !- |
---|
| 378 | zlev_vec(:,:) = lev(1) |
---|
| 379 | !- |
---|
| 380 | !- In most forcing datas obtained from observations, the wind will be available |
---|
| 381 | !- another height than temperature and humidity. We will thus need to lower the |
---|
| 382 | !- wind to take this effect into account. |
---|
| 383 | !- |
---|
| 384 | lower_wind=.FALSE. |
---|
| 385 | IF ( ABS(lev(1) - levuv(1)) > EPSILON(lev) ) THEN |
---|
| 386 | lower_wind=.TRUE. |
---|
| 387 | zlevuv_vec(:,:) = levuv(1) |
---|
| 388 | ELSE |
---|
| 389 | lower_wind=.FALSE. |
---|
| 390 | zlevuv_vec(:,:) = lev(1) |
---|
| 391 | ENDIF |
---|
| 392 | !- |
---|
| 393 | ! prepares kindex table from the information obtained |
---|
| 394 | ! from the forcing data and reads some initial values for |
---|
| 395 | ! temperature, etc. |
---|
| 396 | !- |
---|
| 397 | kindex(1) = 1 |
---|
| 398 | !- |
---|
| 399 | CALL forcing_READ & |
---|
| 400 | & (filename, rest_id, .TRUE., .FALSE., & |
---|
| 401 | & 0, itau_dep, 0, split, nb_spread, netrad_cons, & |
---|
| 402 | & date0, dt_force, iim, jjm, lon, lat, zlev_vec, zlevuv_vec, tm, & |
---|
| 403 | & swdown, precip_rain, precip_snow, tair_obs, & |
---|
| 404 | & u, v, qair_obs, pb, lwdown, for_contfrac, for_neighbours, for_resolution, & |
---|
| 405 | & for_swnet, eair_obs, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, for_ccanopy, & |
---|
| 406 | & kindex, nbindex, force_id) |
---|
| 407 | !- |
---|
| 408 | WRITE (numout,*) ">> Number of land points =",nbindex |
---|
| 409 | IF (nbindex == 0) THEN |
---|
| 410 | WRITE(numout,*) "Limits : (W,E / N,S)", limit_west, limit_east, & |
---|
| 411 | & limit_north, limit_south |
---|
| 412 | CALL ipslerr ( 3, 'dim2_driver','number of land points error.', & |
---|
| 413 | & ' is zero !','stop driver') |
---|
| 414 | ENDIF |
---|
| 415 | |
---|
| 416 | jlandindex = (((kindex(1:nbindex)-1)/iim) + 1) |
---|
| 417 | ilandindex = (kindex(1:nbindex) - (jlandindex(1:nbindex)-1)*iim) |
---|
| 418 | IF (longprint) THEN |
---|
| 419 | WRITE(numout,*) "kindex of land points : ", kindex(1:nbindex) |
---|
| 420 | WRITE(numout,*) "index i of land points : ", ilandindex |
---|
| 421 | WRITE(numout,*) "index j of land points : ", jlandindex |
---|
| 422 | ENDIF |
---|
| 423 | |
---|
| 424 | IF (is_watchout) THEN |
---|
| 425 | lower_wind = .FALSE. |
---|
| 426 | ENDIF |
---|
| 427 | IF (debug) WRITE(numout,*) "lower wind : ",lower_wind |
---|
| 428 | |
---|
| 429 | !- |
---|
| 430 | im = iim; jm = jjm; lm = llm; |
---|
| 431 | IF ( (iim > 1).AND.(jjm > 1) ) THEN |
---|
| 432 | jtest = INT((kindex(INT(nbindex/2))-1)/iim)+1 |
---|
| 433 | itest = MAX( 1, kindex(INT(nbindex/2))-(jtest-1)*iim ) |
---|
| 434 | ELSE |
---|
| 435 | jtest = 1 |
---|
| 436 | itest = 1 |
---|
| 437 | ENDIF |
---|
| 438 | IF (debug) WRITE(numout,*) "test point in dim2_driver : ",itest,jtest |
---|
| 439 | !- |
---|
| 440 | IF ((im /= iim) .AND. (jm /= jjm) .AND. (lm /= llm)) THEN |
---|
| 441 | WRITE(numout,*) ' dimensions are not good. Verify FILE :' |
---|
| 442 | WRITE(numout,*) ' filename = ',filename |
---|
| 443 | WRITE(numout,*) ' im, jm, lm lus = ', im, jm, lm |
---|
| 444 | WRITE(numout,*) ' iim, jjm, llm demandes = ', iim, jjm, llm |
---|
| 445 | STOP 'dim2_driver' |
---|
| 446 | ENDIF |
---|
| 447 | !===================================================================== |
---|
| 448 | !- 1.5 Configures the time-steps and other parameters |
---|
| 449 | !- of the run. |
---|
| 450 | !===================================================================== |
---|
| 451 | !- |
---|
| 452 | ! If the time steping of the restart is different from the one |
---|
| 453 | ! of the forcing we need to convert the itau_dep into the |
---|
| 454 | ! chronology of the forcing. This ensures that the forcing will |
---|
| 455 | ! start at the date of the restart file. Obviously the itau_fin |
---|
| 456 | ! needs to be shifted as well ! |
---|
| 457 | !- |
---|
| 458 | IF ( (dt_rest /= dt_force).AND.(itau_dep > 1) ) THEN |
---|
| 459 | itau_dep = NINT((itau_dep*dt_rest )/dt_force) |
---|
| 460 | itau_fin = itau_dep+tm-1 |
---|
| 461 | if (debug) WRITE(numout,*) & |
---|
| 462 | & 'The time steping of the restart is different from the one ',& |
---|
| 463 | & 'of the forcing we need to convert the itau_dep into the ',& |
---|
| 464 | & 'chronology of the forcing. This ensures that the forcing will ',& |
---|
| 465 | & 'start at the date of the restart file. Obviously the itau_fin ',& |
---|
| 466 | & 'needs to be shifted as well : itau_dep, itau_fin ', & |
---|
| 467 | & itau_dep, itau_fin |
---|
| 468 | ENDIF |
---|
| 469 | !- |
---|
| 470 | ! Same things if the starting dates are not the same. |
---|
| 471 | ! Everything should look as if we had only one forcing file ! |
---|
| 472 | !- |
---|
| 473 | IF (date0 /= date0_rest) THEN |
---|
| 474 | WRITE(numout,*) 'date0_rest , date0 : ',date0_rest , date0 |
---|
| 475 | for_offset = NINT((date0_rest-date0)*one_day/dt_force) |
---|
| 476 | date0 = date0_rest |
---|
| 477 | ELSE |
---|
| 478 | for_offset = 0 |
---|
| 479 | ENDIF |
---|
| 480 | WRITE(numout,*) 'OFFSET FOR THE data read :', for_offset |
---|
| 481 | |
---|
| 482 | CALL ioconf_startdate(date0) |
---|
| 483 | !- |
---|
| 484 | !Config Key = TIME_SKIP |
---|
| 485 | !Config Desc = Time in the forcing file at which the model |
---|
| 486 | !Config is started. |
---|
| 487 | !Config Def = 0 |
---|
| 488 | !Config Help = This time give the point in time at which the model |
---|
| 489 | !Config should be started. If exists, the date of the restart file is use. |
---|
| 490 | !Config The FORMAT of this date can be either of the following : |
---|
| 491 | !Config n : time step n within the forcing file |
---|
| 492 | !Config nS : n seconds after the first time-step in the file |
---|
| 493 | !Config nD : n days after the first time-step |
---|
| 494 | !Config nM : n month after the first time-step (year of 365 days) |
---|
| 495 | !Config nY : n years after the first time-step (year of 365 days) |
---|
| 496 | !Config Or combinations : |
---|
| 497 | !Config nYmM: n years and m month |
---|
| 498 | !- |
---|
| 499 | itau_skip = 0 |
---|
| 500 | WRITE(time_str,'(I10)') itau_skip |
---|
| 501 | CALL getin_p('TIME_SKIP', time_str) |
---|
| 502 | !- |
---|
| 503 | ! Transform into itau |
---|
| 504 | !- |
---|
| 505 | CALL tlen2itau (time_str, dt_force, date0, itau_skip) |
---|
| 506 | !- |
---|
| 507 | itau_dep = itau_dep+itau_skip |
---|
| 508 | !- |
---|
| 509 | ! We need to select the right position of the splited time steps. |
---|
| 510 | !- |
---|
| 511 | istp = itau_dep*split+1 |
---|
| 512 | IF (MOD(istp-1,split) /= 0) THEN |
---|
| 513 | split_start = MOD(istp-1,split)+1 |
---|
| 514 | ELSE |
---|
| 515 | split_start = 1 |
---|
| 516 | ENDIF |
---|
| 517 | istp_old = itau_dep_rest |
---|
| 518 | !!$ it_force = MAX(MOD(itau_dep,INT(one_day*one_year/dt_force)),1)+itau_skip |
---|
| 519 | !!$ WRITE(numout,*) 'itau_dep, it_force :', itau_dep, it_force |
---|
| 520 | !- |
---|
| 521 | itau_len = itau_fin-itau_dep |
---|
| 522 | !- |
---|
| 523 | !Config Key = TIME_LENGTH |
---|
| 524 | !Config Desc = Length of the integration in time. |
---|
| 525 | !Config Def = DEF |
---|
| 526 | !Config Help = Length of integration. By default the entire length |
---|
| 527 | !Config of the forcing is used. The FORMAT of this date can |
---|
| 528 | !Config be either of the following : |
---|
| 529 | !Config n : time step n within the forcing file |
---|
| 530 | !Config nS : n seconds after the first time-step in the file |
---|
| 531 | !Config nD : n days after the first time-step |
---|
| 532 | !Config nM : n month after the first time-step (year of 365 days) |
---|
| 533 | !Config nY : n years after the first time-step (year of 365 days) |
---|
| 534 | !Config Or combinations : |
---|
| 535 | !Config nYmM: n years and m month |
---|
| 536 | !- |
---|
| 537 | WRITE(time_str,'(I10)') itau_len |
---|
| 538 | CALL getin_p('TIME_LENGTH', time_str) |
---|
| 539 | !- |
---|
| 540 | ! Transform into itau |
---|
| 541 | !- |
---|
| 542 | CALL tlen2itau (time_str, dt_force, date0, itau_len) |
---|
| 543 | !- |
---|
| 544 | itau_fin = itau_dep+itau_len |
---|
| 545 | !- |
---|
| 546 | WRITE(numout,*) '>> Time origine in the forcing file :', date0 |
---|
| 547 | WRITE(numout,*) '>> Time origine in the restart file :', date0_rest |
---|
| 548 | WRITE(numout,*) '>> Simulate starts at forcing time-step : ', itau_dep |
---|
| 549 | WRITE(numout,*) '>> The splited time-steps start at (Sets the ' |
---|
| 550 | WRITE(numout,*) '>> chronology for the history and restart files):',istp |
---|
| 551 | WRITE(numout,*) '>> The time spliting starts at :', split_start |
---|
| 552 | WRITE(numout,*) '>> Simulation ends at forcing time-step: ', itau_fin |
---|
| 553 | WRITE(numout,*) '>> Length of the simulation is thus :', itau_len |
---|
| 554 | WRITE(numout,*) '>> Length of the forcing data is in time-steps : ', tm |
---|
| 555 | IF (tm < itau_len) THEN |
---|
| 556 | CALL ipslerr ( 2, 'dim2_driver','Length of the simulation is greater than.', & |
---|
| 557 | & ' Length of the forcing data is in time-steps','verify TIME_LENGTH parameter.') |
---|
| 558 | ENDIF |
---|
| 559 | WRITE(numout,*) '>> Time steps : true, forcing and restart : ', dt,dt_force,dt_rest |
---|
| 560 | |
---|
| 561 | !===================================================================== |
---|
| 562 | !- 2.0 This section is going to define the details by which |
---|
| 563 | !- the input data is going to be used to force the |
---|
| 564 | !- land-surface scheme. The tasks are the following : |
---|
| 565 | !- - Is the coupling going to be explicit or implicit |
---|
| 566 | !- - Type of interpolation to be used. |
---|
| 567 | !- - At which height are the atmospheric forcings going to be used ? |
---|
| 568 | !- - Is a relaxation method going to be used on the forcing |
---|
| 569 | !- - Does net radiation in the interpolated data need to be conserved |
---|
| 570 | !- - How do we distribute precipitation. |
---|
| 571 | !===================================================================== |
---|
| 572 | !Config Key = RELAXATION |
---|
| 573 | !Config Desc = method of forcing |
---|
| 574 | !Config Def = n |
---|
| 575 | !Config Help = A method is proposed by which the first atmospheric |
---|
| 576 | !Config level is not directly forced by observations but |
---|
| 577 | !Config relaxed with a time constant towards observations. |
---|
| 578 | !Config For the moment the methods tends to smooth too much |
---|
| 579 | !Config the diurnal cycle and introduces a time shift. |
---|
| 580 | !Config A more sophisticated method is needed. |
---|
| 581 | !- |
---|
| 582 | relaxation = .FALSE. |
---|
| 583 | CALL getin_p('RELAXATION', relaxation) |
---|
| 584 | IF ( relaxation ) THEN |
---|
| 585 | WRITE(numout,*) 'dim2_driver : The relaxation option is temporarily disabled as it does not' |
---|
| 586 | WRITE(numout,*) ' produce energy conservation in ORCHIDEE. If you intend to use it' |
---|
| 587 | WRITE(numout,*) ' you should validate it.' |
---|
| 588 | STOP 'dim2_driver' |
---|
| 589 | !- |
---|
| 590 | !Config Key = RELAX_A |
---|
| 591 | !Config Desc = Time constant of the relaxation layer |
---|
| 592 | !Config If = RELAXATION |
---|
| 593 | !Config Def = 1.0 |
---|
| 594 | !Config Help = The time constant associated to the atmospheric |
---|
| 595 | !Config conditions which are going to be computed |
---|
| 596 | !Config in the relaxed layer. To avoid too much |
---|
| 597 | !Config damping the value should be larger than 1000. |
---|
| 598 | !- |
---|
| 599 | alpha = 1000.0 |
---|
| 600 | CALL getin_p('RELAX_A', alpha) |
---|
| 601 | ENDIF |
---|
| 602 | !- |
---|
| 603 | !Config Key = NO_INTER||INTER_LIN |
---|
| 604 | !Config Desc = Interpolation or not IF split is larger than 1 |
---|
| 605 | !Config Def = NO_INTER |
---|
| 606 | !Config Help = Choose IF you wish to interpolate linearly or not. |
---|
| 607 | !- |
---|
| 608 | no_inter = .TRUE. |
---|
| 609 | inter_lin = .FALSE. |
---|
| 610 | CALL getin_p('NO_INTER', no_inter) |
---|
| 611 | CALL getin_p('INTER_LIN', inter_lin) |
---|
| 612 | !- |
---|
| 613 | IF (inter_lin) THEN |
---|
| 614 | no_inter = .FALSE. |
---|
| 615 | ELSE |
---|
| 616 | no_inter = .TRUE. |
---|
| 617 | ENDIF |
---|
| 618 | !- |
---|
| 619 | IF (split == 1) THEN |
---|
| 620 | no_inter = .TRUE. |
---|
| 621 | inter_lin = .FALSE. |
---|
| 622 | ENDIF |
---|
| 623 | !- |
---|
| 624 | !Config Key = SPRED_PREC |
---|
| 625 | !Config Desc = Spread the precipitation. |
---|
| 626 | !Config Def = 1 |
---|
| 627 | !Config Help = Spread the precipitation over n steps of the splited forcing |
---|
| 628 | !Config time step. This ONLY applied if the forcing time step has been splited. |
---|
| 629 | !Config If the value indicated is greater than SPLIT_DT, SPLIT_DT is used for it. |
---|
| 630 | !- |
---|
| 631 | nb_spread = 1 |
---|
| 632 | CALL getin_p('SPRED_PREC', nb_spread) |
---|
| 633 | IF (nb_spread > split) THEN |
---|
| 634 | WRITE(numout,*) 'WARNING : nb_spread is too large it will be ' |
---|
| 635 | WRITE(numout,*) ' set to the value of split' |
---|
| 636 | nb_spread = split |
---|
| 637 | ELSE IF (split == 1) THEN |
---|
| 638 | nb_spread = 1 |
---|
| 639 | ENDIF |
---|
| 640 | !- |
---|
| 641 | IF (inter_lin) THEN |
---|
| 642 | !Config Key = NETRAD_CONS |
---|
| 643 | !Config Desc = Conserve net radiation in the forcing |
---|
| 644 | !Config Def = y |
---|
| 645 | !Config If = INTER_LIN |
---|
| 646 | !Config Help = When the interpolation is used the net radiation |
---|
| 647 | !Config provided by the forcing is not conserved anymore. |
---|
| 648 | !Config This should be avoided and thus this option should |
---|
| 649 | !Config be TRUE (y). |
---|
| 650 | !Config This option is not used for short-wave if the |
---|
| 651 | !Config time-step of the forcing is longer than an hour. |
---|
| 652 | !Config It does not make sense to try and reconstruct |
---|
| 653 | !Config a diurnal cycle and at the same time conserve the |
---|
| 654 | !Config incoming solar radiation. |
---|
| 655 | !- |
---|
| 656 | netrad_cons = .TRUE. |
---|
| 657 | CALL getin_p('NETRAD_CONS', netrad_cons) |
---|
| 658 | ELSE |
---|
| 659 | netrad_cons = .FALSE. |
---|
| 660 | ENDIF |
---|
| 661 | !===================================================================== |
---|
| 662 | !- 3.0 Finaly we can prepare all the variables before launching |
---|
| 663 | !- the simulation ! |
---|
| 664 | !===================================================================== |
---|
| 665 | ! Initialize LOGICAL and the length of the integration |
---|
| 666 | !- |
---|
| 667 | last_CALL = .FALSE. |
---|
| 668 | first_CALL = .TRUE. |
---|
| 669 | !- |
---|
| 670 | temp_sol_NEW(:,:) = tp_00 |
---|
| 671 | !- |
---|
| 672 | !Config Key = ATM_CO2 |
---|
| 673 | !Config Desc = Value for atm CO2 |
---|
| 674 | !Config Def = 350. |
---|
| 675 | !Config Help = Value to prescribe the atm CO2. |
---|
| 676 | !Config For pre-industrial simulations, the value is 286.2 . |
---|
| 677 | !Config 348. for 1990 year. |
---|
| 678 | !- |
---|
| 679 | atmco2=350. |
---|
| 680 | CALL getin_p('ATM_CO2',atmco2) |
---|
| 681 | for_ccanopy(:,:)=atmco2 |
---|
| 682 | !- |
---|
| 683 | ! Preparing for the implicit scheme. |
---|
| 684 | ! This means loading the prognostic variables from the restart file. |
---|
| 685 | !- |
---|
| 686 | IF (is_root_prc) & |
---|
| 687 | ALLOCATE(fluxsens_g(iim_g,jjm_g)) |
---|
| 688 | IF (is_root_prc) THEN |
---|
| 689 | var_name= 'fluxsens' |
---|
| 690 | CALL restget & |
---|
| 691 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., fluxsens_g) |
---|
| 692 | IF (ALL(fluxsens_g(:,:) == val_exp)) THEN |
---|
| 693 | fluxsens_g(:,:) = zero |
---|
| 694 | ENDIF |
---|
| 695 | ENDIF |
---|
| 696 | CALL scatter2D(fluxsens_g,fluxsens) |
---|
| 697 | IF (is_root_prc) & |
---|
| 698 | DEALLOCATE(fluxsens_g) |
---|
| 699 | !- |
---|
| 700 | IF (is_root_prc) & |
---|
| 701 | ALLOCATE(vevapp_g(iim_g,jjm_g)) |
---|
| 702 | IF (is_root_prc) THEN |
---|
| 703 | var_name= 'vevapp' |
---|
| 704 | CALL restget & |
---|
| 705 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., vevapp_g) |
---|
| 706 | IF (ALL(vevapp_g(:,:) == val_exp)) THEN |
---|
| 707 | vevapp(:,:) = 0. |
---|
| 708 | ENDIF |
---|
| 709 | ENDIF |
---|
| 710 | CALL scatter2D(vevapp_g,vevapp) |
---|
| 711 | IF (is_root_prc) & |
---|
| 712 | DEALLOCATE(vevapp_g) |
---|
| 713 | !- |
---|
| 714 | IF (is_root_prc) & |
---|
| 715 | ALLOCATE(old_zlev_g(iim_g,jjm_g)) |
---|
| 716 | IF (is_root_prc) THEN |
---|
| 717 | var_name= 'zlev_old' |
---|
| 718 | CALL restget & |
---|
| 719 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., old_zlev_g) |
---|
| 720 | IF (ALL(old_zlev_g(:,:) == val_exp)) THEN |
---|
| 721 | Flag=.TRUE. |
---|
| 722 | ELSE |
---|
| 723 | Flag=.FALSE. |
---|
| 724 | ENDIF |
---|
| 725 | ENDIF |
---|
| 726 | CALL scatter2D(old_zlev_g,old_zlev) |
---|
| 727 | IF (is_root_prc) & |
---|
| 728 | DEALLOCATE(old_zlev_g) |
---|
| 729 | CALL bcast(Flag) |
---|
| 730 | IF ( Flag ) old_zlev(:,:)=zlev_vec(:,:) |
---|
| 731 | !- |
---|
| 732 | IF (is_root_prc) & |
---|
| 733 | ALLOCATE(old_qair_g(iim_g,jjm_g)) |
---|
| 734 | IF (is_root_prc) THEN |
---|
| 735 | var_name= 'qair_old' |
---|
| 736 | CALL restget & |
---|
| 737 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., old_qair_g) |
---|
| 738 | IF (ALL(old_qair_g(:,:) == val_exp)) THEN |
---|
| 739 | Flag=.TRUE. |
---|
| 740 | ELSE |
---|
| 741 | Flag=.FALSE. |
---|
| 742 | ENDIF |
---|
| 743 | ENDIF |
---|
| 744 | CALL scatter2D(old_qair_g,old_qair) |
---|
| 745 | IF (is_root_prc) & |
---|
| 746 | DEALLOCATE(old_qair_g) |
---|
| 747 | CALL bcast(Flag) |
---|
| 748 | IF (Flag) old_qair(:,:) = qair_obs(:,:) |
---|
| 749 | !- |
---|
| 750 | IF (is_root_prc) & |
---|
| 751 | ALLOCATE(old_eair_g(iim_g,jjm_g)) |
---|
| 752 | IF (is_root_prc) THEN |
---|
| 753 | var_name= 'eair_old' |
---|
| 754 | CALL restget & |
---|
| 755 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., old_eair_g) |
---|
| 756 | IF (ALL(old_eair_g(:,:) == val_exp)) THEN |
---|
| 757 | Flag=.TRUE. |
---|
| 758 | ELSE |
---|
| 759 | Flag=.FALSE. |
---|
| 760 | ENDIF |
---|
| 761 | ENDIF |
---|
| 762 | CALL scatter2D(old_eair_g,old_eair) |
---|
| 763 | IF (is_root_prc) & |
---|
| 764 | DEALLOCATE(old_eair_g) |
---|
| 765 | CALL bcast(Flag) |
---|
| 766 | IF (Flag) THEN |
---|
| 767 | DO ik=1,nbindex |
---|
| 768 | i=ilandindex(ik) |
---|
| 769 | j=jlandindex(ik) |
---|
| 770 | old_eair(i,j) = cp_air * tair_obs(i,j) + cte_grav*zlev_vec(i,j) |
---|
| 771 | ENDDO |
---|
| 772 | ENDIF |
---|
| 773 | !- |
---|
| 774 | ! old density is also needed because we do not yet have the right pb |
---|
| 775 | !- |
---|
| 776 | !=> obsolète ??!! (tjrs calculé après forcing_read) |
---|
| 777 | IF (is_root_prc) & |
---|
| 778 | ALLOCATE(for_rau_g(iim_g,jjm_g)) |
---|
| 779 | IF (is_root_prc) THEN |
---|
| 780 | var_name= 'rau_old' |
---|
| 781 | CALL restget & |
---|
| 782 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., for_rau_g) |
---|
| 783 | IF (ALL(for_rau_g(:,:) == val_exp)) THEN |
---|
| 784 | Flag=.TRUE. |
---|
| 785 | ELSE |
---|
| 786 | Flag=.FALSE. |
---|
| 787 | ENDIF |
---|
| 788 | ENDIF |
---|
| 789 | CALL scatter2D(for_rau_g,for_rau) |
---|
| 790 | IF (is_root_prc) & |
---|
| 791 | DEALLOCATE(for_rau_g) |
---|
| 792 | CALL bcast(Flag) |
---|
| 793 | IF (Flag) THEN |
---|
| 794 | DO ik=1,nbindex |
---|
| 795 | i=ilandindex(ik) |
---|
| 796 | j=jlandindex(ik) |
---|
| 797 | for_rau(i,j) = pb(i,j)/(cte_molr*(tair_obs(i,j))) |
---|
| 798 | ENDDO |
---|
| 799 | ENDIF |
---|
| 800 | !- |
---|
| 801 | ! For this variable the restart is extracted by SECHIBA |
---|
| 802 | !- |
---|
| 803 | temp_sol_NEW(:,:) = tair_obs(:,:) |
---|
| 804 | !- |
---|
| 805 | if (.NOT. is_watchout) THEN |
---|
| 806 | !- |
---|
| 807 | ! This does not yield a correct restart in the case of relaxation |
---|
| 808 | !- |
---|
| 809 | IF (is_root_prc) & |
---|
| 810 | ALLOCATE(petAcoef_g(iim_g,jjm_g)) |
---|
| 811 | IF (is_root_prc) THEN |
---|
| 812 | var_name= 'petAcoef' |
---|
| 813 | CALL restget & |
---|
| 814 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., petAcoef_g) |
---|
| 815 | IF (ALL(petAcoef_g(:,:) == val_exp)) THEN |
---|
| 816 | Flag=.TRUE. |
---|
| 817 | ELSE |
---|
| 818 | Flag=.FALSE. |
---|
| 819 | ENDIF |
---|
| 820 | ENDIF |
---|
| 821 | CALL scatter2D(petAcoef_g,petAcoef) |
---|
| 822 | IF (is_root_prc) & |
---|
| 823 | DEALLOCATE(petAcoef_g) |
---|
| 824 | CALL bcast(Flag) |
---|
| 825 | IF (Flag) petAcoef(:,:) = zero |
---|
| 826 | !-- |
---|
| 827 | IF (is_root_prc) & |
---|
| 828 | ALLOCATE(petBcoef_g(iim_g,jjm_g)) |
---|
| 829 | IF (is_root_prc) THEN |
---|
| 830 | var_name= 'petBcoef' |
---|
| 831 | CALL restget & |
---|
| 832 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., petBcoef_g) |
---|
| 833 | IF (ALL(petBcoef_g(:,:) == val_exp)) THEN |
---|
| 834 | Flag=.TRUE. |
---|
| 835 | ELSE |
---|
| 836 | Flag=.FALSE. |
---|
| 837 | ENDIF |
---|
| 838 | ENDIF |
---|
| 839 | CALL scatter2D(petBcoef_g,petBcoef) |
---|
| 840 | IF (is_root_prc) & |
---|
| 841 | DEALLOCATE(petBcoef_g) |
---|
| 842 | CALL bcast(Flag) |
---|
| 843 | IF (Flag) petBcoef(:,:) = old_eair(:,:) |
---|
| 844 | !-- |
---|
| 845 | IF (is_root_prc) & |
---|
| 846 | ALLOCATE(peqAcoef_g(iim_g,jjm_g)) |
---|
| 847 | IF (is_root_prc) THEN |
---|
| 848 | var_name= 'peqAcoef' |
---|
| 849 | CALL restget & |
---|
| 850 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., peqAcoef_g) |
---|
| 851 | IF (ALL(peqAcoef_g(:,:) == val_exp)) THEN |
---|
| 852 | Flag=.TRUE. |
---|
| 853 | ELSE |
---|
| 854 | Flag=.FALSE. |
---|
| 855 | ENDIF |
---|
| 856 | ENDIF |
---|
| 857 | CALL scatter2D(peqAcoef_g,peqAcoef) |
---|
| 858 | IF (is_root_prc) & |
---|
| 859 | DEALLOCATE(peqAcoef_g) |
---|
| 860 | CALL bcast(Flag) |
---|
| 861 | IF (Flag) peqAcoef(:,:) = zero |
---|
| 862 | !-- |
---|
| 863 | IF (is_root_prc) & |
---|
| 864 | ALLOCATE(peqBcoef_g(iim_g,jjm_g)) |
---|
| 865 | IF (is_root_prc) THEN |
---|
| 866 | var_name= 'peqBcoef' |
---|
| 867 | CALL restget & |
---|
| 868 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., peqBcoef_g) |
---|
| 869 | IF (ALL(peqBcoef_g(:,:) == val_exp)) THEN |
---|
| 870 | Flag=.TRUE. |
---|
| 871 | ELSE |
---|
| 872 | Flag=.FALSE. |
---|
| 873 | ENDIF |
---|
| 874 | ENDIF |
---|
| 875 | CALL scatter2D(peqBcoef_g,peqBcoef) |
---|
| 876 | IF (is_root_prc) & |
---|
| 877 | DEALLOCATE(peqBcoef_g) |
---|
| 878 | CALL bcast(Flag) |
---|
| 879 | IF (Flag) peqBcoef(:,:) = old_qair(:,:) |
---|
| 880 | ENDIF |
---|
| 881 | !- |
---|
| 882 | ! And other variables which need initial variables. These variables |
---|
| 883 | ! will get properly initialized by ORCHIDEE when it is called for |
---|
| 884 | ! the first time. |
---|
| 885 | !- |
---|
| 886 | albedo(:,:,:) = 0.13 |
---|
| 887 | emis(:,:) = 1.0 |
---|
| 888 | z0(:,:) = 0.1 |
---|
| 889 | !-- |
---|
| 890 | !===================================================================== |
---|
| 891 | !- 4.0 START THE TIME LOOP |
---|
| 892 | !===================================================================== |
---|
| 893 | |
---|
| 894 | CALL init_timer |
---|
| 895 | |
---|
| 896 | it = itau_dep+1 |
---|
| 897 | DO WHILE ( it <= itau_fin ) |
---|
| 898 | !---- |
---|
| 899 | it_force = it+for_offset |
---|
| 900 | IF (it_force < 0) THEN |
---|
| 901 | WRITE(numout,*) 'The day is not in the forcing file :', & |
---|
| 902 | & it_force, it, for_offset |
---|
| 903 | STOP 'dim2_driver' |
---|
| 904 | ENDIF |
---|
| 905 | !!$ IF (it_force > itau_dep+tm) THEN |
---|
| 906 | !!$ WRITE(numout,*) 'ERROR : more time-steps than data' |
---|
| 907 | !!$ WRITE(numout,*) 'it_force : ', it_force, ' itau_dep+tm : ', itau_dep+tm |
---|
| 908 | !!$ STOP 'dim2_driver' |
---|
| 909 | !!$ ENDIF |
---|
| 910 | !--- |
---|
| 911 | is=split_start |
---|
| 912 | DO WHILE ( is <= split ) |
---|
| 913 | !----- |
---|
| 914 | julian = itau2date(istp, date0, dt) |
---|
| 915 | CALL ju2ymds(julian, yy, mm, dd, ss) |
---|
| 916 | IF (debug) THEN |
---|
| 917 | WRITE(numout,*) "==============================================================" |
---|
| 918 | WRITE(numout,"('New iteration at date : ',I4,'-',I2.2,'-',I2.2,':',F8.4)") & |
---|
| 919 | & yy,mm,dd,ss/3600. |
---|
| 920 | #ifdef CPP_PARA |
---|
| 921 | IF (is_root_prc) THEN |
---|
| 922 | WRITE(*,*) "==============================================================" |
---|
| 923 | WRITE(*,"('New iteration at date : ',I4,'-',I2.2,'-',I2.2,':',F8.4)") & |
---|
| 924 | & yy,mm,dd,ss/3600. |
---|
| 925 | ENDIF |
---|
| 926 | #endif |
---|
| 927 | ENDIF |
---|
| 928 | !----- |
---|
| 929 | IF ( (it == itau_fin).AND.(is == split) ) THEN |
---|
| 930 | last_CALL = .TRUE. |
---|
| 931 | ENDIF |
---|
| 932 | !----- |
---|
| 933 | IF (debug) WRITE(numout,*) 'Into forcing_read' |
---|
| 934 | !----- |
---|
| 935 | CALL forcing_READ & |
---|
| 936 | & (filename, rest_id, .FALSE., last_call, & |
---|
| 937 | & it_force, istp, is, split, nb_spread, netrad_cons, & |
---|
| 938 | & date0, dt_force, iim, jjm, lon, lat, zlev_vec, zlevuv_vec, tm, & |
---|
| 939 | & swdown, precip_rain, precip_snow, tair_obs, & |
---|
| 940 | & u, v, qair_obs, pb, for_lwdown, for_contfrac, for_neighbours, for_resolution, & |
---|
| 941 | & for_swnet, eair_obs, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, for_ccanopy, & |
---|
| 942 | & kindex, nbindex, force_id) |
---|
| 943 | |
---|
| 944 | !----- |
---|
| 945 | !---- SECHIBA expects surface pressure in hPa |
---|
| 946 | !----- |
---|
| 947 | for_psurf(:,:) = pb(:,:)/100. |
---|
| 948 | |
---|
| 949 | IF (longprint) THEN |
---|
| 950 | WRITE(numout,*) "dim2_driver 0 ",it_force |
---|
| 951 | WRITE(numout,*) ">> Index of land points =",kindex |
---|
| 952 | WRITE(numout,*) "Lowest level wind speed North = ", & |
---|
| 953 | & (/ ( u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 954 | WRITE(numout,*) "Lowest level wind speed East = ", & |
---|
| 955 | & (/ ( v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 956 | WRITE(numout,*) "z0 ; Surface roughness = ", & |
---|
| 957 | & (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 958 | WRITE(numout,*) "Height of first layer = ", & |
---|
| 959 | & (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 960 | WRITE(numout,*) "Lowest level specific humidity = ", & |
---|
| 961 | & (/ ( qair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 962 | WRITE(numout,*) "Rain precipitation = ", & |
---|
| 963 | & (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /) |
---|
| 964 | WRITE(numout,*) "Snow precipitation = ", & |
---|
| 965 | & (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /) |
---|
| 966 | WRITE(numout,*) "Down-welling long-wave flux = ", & |
---|
| 967 | & (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 968 | WRITE(numout,*) "Net surface short-wave flux = ", & |
---|
| 969 | & (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 970 | WRITE(numout,*) "Downwelling surface short-wave flux = ", & |
---|
| 971 | & (/ ( swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 972 | WRITE(numout,*) "Air temperature in Kelvin = ", & |
---|
| 973 | & (/ ( tair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 974 | WRITE(numout,*) "Air potential energy = ", & |
---|
| 975 | & (/ ( eair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 976 | WRITE(numout,*) "CO2 concentration in the canopy = ", & |
---|
| 977 | & (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 978 | WRITE(numout,*) "Coeficients A from the PBL resolution = ", & |
---|
| 979 | & (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 980 | WRITE(numout,*) "One for T and another for q = ", & |
---|
| 981 | & (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 982 | WRITE(numout,*) "Coeficients B from the PBL resolution = ", & |
---|
| 983 | & (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 984 | WRITE(numout,*) "One for T and another for q = ", & |
---|
| 985 | & (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 986 | WRITE(numout,*) "Cdrag = ", & |
---|
| 987 | & (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 988 | WRITE(numout,*) "CO2 concentration in the canopy = ", & |
---|
| 989 | & (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 990 | WRITE(numout,*) "Lowest level pressure = ", & |
---|
| 991 | & (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 992 | WRITE(numout,*) "Geographical coordinates lon = ", & |
---|
| 993 | & (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 994 | WRITE(numout,*) "Geographical coordinates lat = ", & |
---|
| 995 | & (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 996 | WRITE(numout,*) "Fraction of continent in the grid = ", & |
---|
| 997 | & (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 998 | ENDIF |
---|
| 999 | !----- |
---|
| 1000 | !---- Prepare : tmp_qair, tmp_eair, tmp_tair, tmp_pb |
---|
| 1001 | !---- and : for_u, for_v, for_lwdown, for_swnet, for_swdown |
---|
| 1002 | !---- All the work is done in forcing_read |
---|
| 1003 | !---- We do not need the no_inter options as we will |
---|
| 1004 | !---- allways interpolate |
---|
| 1005 | !----- |
---|
| 1006 | IF (debug) WRITE(numout,*) 'Prepare the atmospheric forcing' |
---|
| 1007 | !----- |
---|
| 1008 | IF (.NOT. is_watchout) THEN |
---|
| 1009 | DO ik=1,nbindex |
---|
| 1010 | i=ilandindex(ik) |
---|
| 1011 | j=jlandindex(ik) |
---|
| 1012 | eair_obs(i,j) = cp_air*tair_obs(i,j)+cte_grav*zlev_vec(i,j) |
---|
| 1013 | for_swnet(i,j) = (1.-(albedo(i,j,1)+albedo(i,j,2))/2.)*swdown(i,j) |
---|
| 1014 | ENDDO |
---|
| 1015 | ENDIF |
---|
| 1016 | for_swdown(:,:) = swdown(:,:) |
---|
| 1017 | !----- |
---|
| 1018 | !---- Computing the buffer zone ! |
---|
| 1019 | !----- |
---|
| 1020 | IF (relaxation) THEN |
---|
| 1021 | DO ik=1,nbindex |
---|
| 1022 | i=ilandindex(ik) |
---|
| 1023 | j=jlandindex(ik) |
---|
| 1024 | for_qair(i,j) = peqAcoef(i,j)*(-1.) * vevapp(i,j)*dt+peqBcoef(i,j) |
---|
| 1025 | !------- |
---|
| 1026 | for_eair(i,j) = petAcoef(i,j)*(-1.) * fluxsens(i,j)+petBcoef(i,j) |
---|
| 1027 | !------- |
---|
| 1028 | ENDDO |
---|
| 1029 | DO ik=1,nbindex |
---|
| 1030 | i=ilandindex(ik) |
---|
| 1031 | j=jlandindex(ik) |
---|
| 1032 | for_tair(i,j) = (for_eair(i,j) - cte_grav*zlev_vec(i,j))/cp_air |
---|
| 1033 | !------- |
---|
| 1034 | !!$ if (.NOT. is_watchout) & |
---|
| 1035 | !!$ epot_sol(:,:) = cp_air*temp_sol_NEW(:,:) |
---|
| 1036 | !------- |
---|
| 1037 | ENDDO |
---|
| 1038 | DO ik=1,nbindex |
---|
| 1039 | i=ilandindex(ik) |
---|
| 1040 | j=jlandindex(ik) |
---|
| 1041 | for_rau(i,j) = pb(i,j) / (cte_molr*for_tair(i,j)) |
---|
| 1042 | !------- |
---|
| 1043 | relax(i,j) = for_rau(i,j)*alpha |
---|
| 1044 | ENDDO |
---|
| 1045 | zlflu = lev(1)/2.0*dt |
---|
| 1046 | DO ik=1,nbindex |
---|
| 1047 | i=ilandindex(ik) |
---|
| 1048 | j=jlandindex(ik) |
---|
| 1049 | peqAcoef(i,j) = 1.0/(zlflu+relax(i,j)) |
---|
| 1050 | peqBcoef(i,j) = (relax(i,j) * qair_obs(i,j)/(zlflu+relax(i,j))) + & |
---|
| 1051 | & for_qair(i,j)/(1.0+relax(i,j)/zlflu) |
---|
| 1052 | ENDDO |
---|
| 1053 | !------- |
---|
| 1054 | ! relax(:,:) = for_rau(:,:)*alpha |
---|
| 1055 | DO ik=1,nbindex |
---|
| 1056 | i=ilandindex(ik) |
---|
| 1057 | j=jlandindex(ik) |
---|
| 1058 | petAcoef(i,j) = 1.0/(zlflu+relax(i,j)) |
---|
| 1059 | petBcoef(i,j) = ( relax(i,j) * eair_obs(i,j) / (zlflu+relax(i,j)) ) & |
---|
| 1060 | & + for_eair(i,j)/(1.0+relax(i,j)/zlflu) |
---|
| 1061 | ENDDO |
---|
| 1062 | ELSE |
---|
| 1063 | for_qair(:,:) = fill_init |
---|
| 1064 | for_eair(:,:) = fill_init |
---|
| 1065 | for_tair(:,:) = fill_init |
---|
| 1066 | DO ik=1,nbindex |
---|
| 1067 | i=ilandindex(ik) |
---|
| 1068 | j=jlandindex(ik) |
---|
| 1069 | for_qair(i,j) = qair_obs(i,j) |
---|
| 1070 | for_eair(i,j) = eair_obs(i,j) |
---|
| 1071 | for_tair(i,j) = tair_obs(i,j) |
---|
| 1072 | ENDDO |
---|
| 1073 | !------- |
---|
| 1074 | !!$ if (.NOT. is_watchout) & |
---|
| 1075 | !!$ epot_sol(:,:) = cp_air*temp_sol_NEW(:,:) |
---|
| 1076 | !------- |
---|
| 1077 | DO ik=1,nbindex |
---|
| 1078 | i=ilandindex(ik) |
---|
| 1079 | j=jlandindex(ik) |
---|
| 1080 | for_rau(i,j) = pb(i,j) / (cte_molr*for_tair(i,j)) |
---|
| 1081 | ENDDO |
---|
| 1082 | !------- |
---|
| 1083 | IF (.NOT. is_watchout) THEN |
---|
| 1084 | petAcoef(:,:) = 0.0 |
---|
| 1085 | peqAcoef(:,:) = 0.0 |
---|
| 1086 | DO ik=1,nbindex |
---|
| 1087 | i=ilandindex(ik) |
---|
| 1088 | j=jlandindex(ik) |
---|
| 1089 | petBcoef(i,j) = eair_obs(i,j) |
---|
| 1090 | peqBcoef(i,j) = qair_obs(i,j) |
---|
| 1091 | ENDDO |
---|
| 1092 | ENDIF |
---|
| 1093 | ENDIF |
---|
| 1094 | !----- |
---|
| 1095 | IF (.NOT. is_watchout) & |
---|
| 1096 | cdrag(:,:) = 0.0 |
---|
| 1097 | for_ccanopy(:,:)=atmco2 |
---|
| 1098 | !----- |
---|
| 1099 | !---- SECHIBA expects wind, temperature and humidity at the same height. |
---|
| 1100 | !---- If this is not the case then we need to correct for that. |
---|
| 1101 | !----- |
---|
| 1102 | IF ( lower_wind .AND. .NOT. is_watchout ) THEN |
---|
| 1103 | DO ik=1,nbindex |
---|
| 1104 | i=ilandindex(ik) |
---|
| 1105 | j=jlandindex(ik) |
---|
| 1106 | for_u(i,j) = u(i,j)*LOG(zlev_vec(i,j)/z0(i,j)) / & |
---|
| 1107 | & LOG(zlevuv_vec(i,j)/z0(i,j)) |
---|
| 1108 | for_v(i,j) = v(i,j)*LOG(zlev_vec(i,j)/z0(i,j)) / & |
---|
| 1109 | & LOG(zlevuv_vec(i,j)/z0(i,j)) |
---|
| 1110 | ENDDO |
---|
| 1111 | ELSE |
---|
| 1112 | DO ik=1,nbindex |
---|
| 1113 | i=ilandindex(ik) |
---|
| 1114 | j=jlandindex(ik) |
---|
| 1115 | for_u(i,j) = u(i,j) |
---|
| 1116 | for_v(i,j) = v(i,j) |
---|
| 1117 | ENDDO |
---|
| 1118 | ENDIF |
---|
| 1119 | !----- |
---|
| 1120 | !---- Prepare the other variables WITH the special CASE |
---|
| 1121 | !---- of splited time steps |
---|
| 1122 | !---- |
---|
| 1123 | !---- PRINT input value for debug |
---|
| 1124 | !----- |
---|
| 1125 | IF (debug) THEN |
---|
| 1126 | WRITE(numout,*) ' >>>>>> time step it_force = ',it_force |
---|
| 1127 | WRITE(numout,*) & |
---|
| 1128 | & ' tair, qair, eair = ', & |
---|
| 1129 | & for_tair(itest,jtest),for_qair(itest,jtest), & |
---|
| 1130 | & for_eair(itest,jtest) |
---|
| 1131 | WRITE(numout,*) & |
---|
| 1132 | & ' OBS : tair, qair, eair = ', & |
---|
| 1133 | & tair_obs(itest,jtest),qair_obs(itest,jtest), & |
---|
| 1134 | & eair_obs(itest,jtest) |
---|
| 1135 | WRITE(numout,*) ' u et v = ',for_u(itest,jtest),for_v(itest,jtest) |
---|
| 1136 | WRITE(numout,*) ' precip rain et snow = ', & |
---|
| 1137 | & precip_rain(itest,jtest),precip_snow(itest,jtest) |
---|
| 1138 | WRITE(numout,*) ' lwdown et swnet = ', & |
---|
| 1139 | & for_lwdown(itest,jtest),for_swnet(itest,jtest) |
---|
| 1140 | WRITE(numout,*) ' petAcoef et peqAcoef = ', & |
---|
| 1141 | & petAcoef(itest,jtest), peqAcoef(itest,jtest) |
---|
| 1142 | WRITE(numout,*) ' petBcoef et peqAcoef = ', & |
---|
| 1143 | & petBcoef(itest,jtest),peqBcoef(itest,jtest) |
---|
| 1144 | WRITE(numout,*) ' zlev = ',zlev_vec(itest,jtest) |
---|
| 1145 | ENDIF |
---|
| 1146 | !----- |
---|
| 1147 | IF (first_CALL) THEN |
---|
| 1148 | |
---|
| 1149 | DO ik=1,nbindex |
---|
| 1150 | i=ilandindex(ik) |
---|
| 1151 | j=jlandindex(ik) |
---|
| 1152 | for_swdown(i,j) = swdown(i,j) |
---|
| 1153 | ENDDO |
---|
| 1154 | IF (longprint) THEN |
---|
| 1155 | WRITE(numout,*) "dim2_driver first_CALL ",it_force |
---|
| 1156 | WRITE(numout,*) ">> Index of land points =",kindex |
---|
| 1157 | WRITE(numout,*) "Lowest level wind speed North = ", & |
---|
| 1158 | & (/ ( for_u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1159 | WRITE(numout,*) "Lowest level wind speed East = ", & |
---|
| 1160 | & (/ ( for_v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1161 | WRITE(numout,*) "z0 ; Surface roughness = ", & |
---|
| 1162 | & (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1163 | WRITE(numout,*) "Height of first layer = ", & |
---|
| 1164 | & (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1165 | WRITE(numout,*) "Lowest level specific humidity = ", & |
---|
| 1166 | & (/ ( for_qair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1167 | WRITE(numout,*) "Rain precipitation = ", & |
---|
| 1168 | & (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /) |
---|
| 1169 | WRITE(numout,*) "Snow precipitation = ", & |
---|
| 1170 | & (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /) |
---|
| 1171 | WRITE(numout,*) "Down-welling long-wave flux = ", & |
---|
| 1172 | & (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1173 | WRITE(numout,*) "Net surface short-wave flux = ", & |
---|
| 1174 | & (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1175 | WRITE(numout,*) "Downwelling surface short-wave flux = ", & |
---|
| 1176 | & (/ ( for_swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1177 | WRITE(numout,*) "Air temperature in Kelvin = ", & |
---|
| 1178 | & (/ ( for_tair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1179 | WRITE(numout,*) "Air potential energy = ", & |
---|
| 1180 | & (/ ( for_eair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1181 | WRITE(numout,*) "CO2 concentration in the canopy = ", & |
---|
| 1182 | & (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1183 | WRITE(numout,*) "Coeficients A from the PBL resolution = ", & |
---|
| 1184 | & (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1185 | WRITE(numout,*) "One for T and another for q = ", & |
---|
| 1186 | & (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1187 | WRITE(numout,*) "Coeficients B from the PBL resolution = ", & |
---|
| 1188 | & (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1189 | WRITE(numout,*) "One for T and another for q = ", & |
---|
| 1190 | & (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1191 | WRITE(numout,*) "Cdrag = ", & |
---|
| 1192 | & (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1193 | WRITE(numout,*) "Lowest level pressure = ", & |
---|
| 1194 | & (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1195 | WRITE(numout,*) "Geographical coordinates lon = ", & |
---|
| 1196 | & (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1197 | WRITE(numout,*) "Geographical coordinates lat = ", & |
---|
| 1198 | & (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1199 | WRITE(numout,*) "Fraction of continent in the grid = ", & |
---|
| 1200 | & (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1201 | ENDIF |
---|
| 1202 | !------- |
---|
| 1203 | !------ CALL sechiba to initialize fields |
---|
| 1204 | !------ and have some initial results: emis, albedo, z0 |
---|
| 1205 | !------- |
---|
| 1206 | CALL intersurf_main & |
---|
| 1207 | & (istp_old, iim, jjm, nbindex, kindex, dt, & |
---|
| 1208 | & first_CALL, .FALSE., lon, lat, for_contfrac, for_neighbours, for_resolution, date0, & |
---|
| 1209 | ! first level conditions |
---|
| 1210 | & zlev_vec, for_u, for_v, & |
---|
| 1211 | & for_qair, for_tair, for_eair, for_ccanopy, & |
---|
| 1212 | ! Variables for the implicit coupling |
---|
| 1213 | & cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, & |
---|
| 1214 | ! Rain, snow, radiation and surface pressure |
---|
| 1215 | & precip_rain, precip_snow, & |
---|
| 1216 | & for_lwdown, for_swnet, for_swdown, for_psurf, & |
---|
| 1217 | ! Output : Fluxes |
---|
| 1218 | & vevapp, fluxsens, fluxlat, coastalflow, riverflow, & |
---|
| 1219 | ! Surface temperatures and surface properties |
---|
| 1220 | & tsol_rad, temp_sol_NEW, qsurf, albedo, emis, z0) |
---|
| 1221 | |
---|
| 1222 | ! |
---|
| 1223 | first_CALL = .FALSE. |
---|
| 1224 | ! |
---|
| 1225 | ! Get Restart values for albedo and z0, |
---|
| 1226 | ! as they modify forcing variables swnet and wind. |
---|
| 1227 | !------- |
---|
| 1228 | ! albedo |
---|
| 1229 | IF (is_root_prc) & |
---|
| 1230 | ALLOCATE(albedo_g(iim_g,jjm_g)) |
---|
| 1231 | ! |
---|
| 1232 | IF (is_root_prc) THEN |
---|
| 1233 | var_name= 'albedo_vis' |
---|
| 1234 | CALL restget & |
---|
| 1235 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., albedo_g) |
---|
| 1236 | IF (ALL(albedo_g(:,:) == val_exp)) THEN |
---|
| 1237 | Flag=.TRUE. |
---|
| 1238 | ELSE |
---|
| 1239 | Flag=.FALSE. |
---|
| 1240 | ENDIF |
---|
| 1241 | ENDIF |
---|
| 1242 | CALL scatter2D(albedo_g,albedo_vis) |
---|
| 1243 | CALL bcast(Flag) |
---|
| 1244 | IF (.NOT. Flag) albedo(:,:,1)=albedo_vis(:,:) |
---|
| 1245 | ! |
---|
| 1246 | IF (is_root_prc) THEN |
---|
| 1247 | var_name= 'albedo_nir' |
---|
| 1248 | CALL restget & |
---|
| 1249 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., albedo_g) |
---|
| 1250 | IF (ALL(albedo_g(:,:) == val_exp)) THEN |
---|
| 1251 | Flag=.TRUE. |
---|
| 1252 | ELSE |
---|
| 1253 | Flag=.FALSE. |
---|
| 1254 | ENDIF |
---|
| 1255 | ENDIF |
---|
| 1256 | CALL scatter2D(albedo_g,albedo_nir) |
---|
| 1257 | CALL bcast(Flag) |
---|
| 1258 | IF (.NOT. Flag) albedo(:,:,2)=albedo_nir(:,:) |
---|
| 1259 | ! |
---|
| 1260 | IF (is_root_prc) & |
---|
| 1261 | DEALLOCATE(albedo_g) |
---|
| 1262 | !-- |
---|
| 1263 | ! z0 |
---|
| 1264 | IF (is_root_prc) & |
---|
| 1265 | ALLOCATE(z0_g(iim_g,jjm_g)) |
---|
| 1266 | IF (is_root_prc) THEN |
---|
| 1267 | var_name= 'z0' |
---|
| 1268 | CALL restget & |
---|
| 1269 | & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., z0_g) |
---|
| 1270 | IF (ALL(z0_g(:,:) == val_exp)) THEN |
---|
| 1271 | Flag=.TRUE. |
---|
| 1272 | ELSE |
---|
| 1273 | Flag=.FALSE. |
---|
| 1274 | ENDIF |
---|
| 1275 | ENDIF |
---|
| 1276 | CALL bcast(Flag) |
---|
| 1277 | IF (.NOT. Flag) CALL scatter2D(z0_g,z0) |
---|
| 1278 | IF (is_root_prc) & |
---|
| 1279 | DEALLOCATE(z0_g) |
---|
| 1280 | !------- |
---|
| 1281 | DO ik=1,nbindex |
---|
| 1282 | i=ilandindex(ik) |
---|
| 1283 | j=jlandindex(ik) |
---|
| 1284 | temp_sol_old(i,j) = temp_sol_NEW(i,j) |
---|
| 1285 | for_swnet(i,j) = (1.- (albedo(i,j,1)+albedo(i,j,2))/2.)*swdown(i,j) |
---|
| 1286 | for_swdown(i,j) = swdown(i,j) |
---|
| 1287 | ENDDO |
---|
| 1288 | ! |
---|
| 1289 | ! MM : z0 have been modified then we must lower the wind again |
---|
| 1290 | !----- |
---|
| 1291 | !---- SECHIBA expects wind, temperature and humidity at the same height. |
---|
| 1292 | !---- If this is not the case then we need to correct for that. |
---|
| 1293 | !----- |
---|
| 1294 | IF ( lower_wind .AND. .NOT. is_watchout ) THEN |
---|
| 1295 | DO ik=1,nbindex |
---|
| 1296 | i=ilandindex(ik) |
---|
| 1297 | j=jlandindex(ik) |
---|
| 1298 | for_u(i,j) = u(i,j) * LOG(zlev_vec(i,j)/z0(i,j)) / & |
---|
| 1299 | & LOG(zlevuv_vec(i,j)/z0(i,j)) |
---|
| 1300 | for_v(i,j) = v(i,j) * LOG(zlev_vec(i,j)/z0(i,j)) / & |
---|
| 1301 | & LOG(zlevuv_vec(i,j)/z0(i,j)) |
---|
| 1302 | ENDDO |
---|
| 1303 | ELSE |
---|
| 1304 | DO ik=1,nbindex |
---|
| 1305 | i=ilandindex(ik) |
---|
| 1306 | j=jlandindex(ik) |
---|
| 1307 | for_u(i,j) = u(i,j) |
---|
| 1308 | for_v(i,j) = v(i,j) |
---|
| 1309 | ENDDO |
---|
| 1310 | ENDIF |
---|
| 1311 | !----- |
---|
| 1312 | !---- PRINT input value after first_CALL for debug |
---|
| 1313 | !----- |
---|
| 1314 | IF (debug) THEN |
---|
| 1315 | WRITE(numout,*) ' >>>>>> after first_CALL = ',first_CALL |
---|
| 1316 | WRITE(numout,*) ' u et v = ',for_u(itest,jtest),for_v(itest,jtest) |
---|
| 1317 | WRITE(numout,*) ' swnet = ', for_swnet(itest,jtest) |
---|
| 1318 | ENDIF |
---|
| 1319 | !------- |
---|
| 1320 | IF (longprint) THEN |
---|
| 1321 | WRITE(numout,*) "dim2_driver first_CALL outputs" |
---|
| 1322 | ! Output : Fluxes |
---|
| 1323 | WRITE(numout,*) "vevapp ; Total of evaporation = ", & |
---|
| 1324 | & (/ ( vevapp(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1325 | WRITE(numout,*) "Sensible heat flux = ", & |
---|
| 1326 | & (/ ( fluxsens(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1327 | WRITE(numout,*) "Latent heat flux = ", & |
---|
| 1328 | & (/ ( fluxlat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1329 | WRITE(numout,*) "coastalflow ; Diffuse flow of water into the ocean (m^3/dt) = ", & |
---|
| 1330 | & (/ ( coastalflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1331 | WRITE(numout,*) "riverflow ; Largest rivers flowing into the ocean (m^3/dt) = ", & |
---|
| 1332 | & (/ ( riverflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1333 | ! Surface temperatures and surface properties |
---|
| 1334 | WRITE(numout,*) "tsol_rad ; Radiative surface temperature = ", & |
---|
| 1335 | & (/ ( tsol_rad(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1336 | WRITE(numout,*) "temp_sol_new ; New soil temperature = ", & |
---|
| 1337 | & (/ ( temp_sol_NEW(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1338 | WRITE(numout,*) "qsurf ; Surface specific humidity = ", & |
---|
| 1339 | & (/ ( qsurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1340 | WRITE(numout,*) "albedoVIS = ", & |
---|
| 1341 | & (/ ( albedo(ilandindex(ik), jlandindex(ik), 1),ik=1,nbindex ) /) |
---|
| 1342 | WRITE(numout,*) "albedoNIR = ", & |
---|
| 1343 | & (/ ( albedo(ilandindex(ik), jlandindex(ik), 2),ik=1,nbindex ) /) |
---|
| 1344 | WRITE(numout,*) "emis ; Emissivity = ", & |
---|
| 1345 | & (/ ( emis(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1346 | WRITE(numout,*) "z0 ; Surface roughness = ", & |
---|
| 1347 | & (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1348 | ENDIF |
---|
| 1349 | !------- |
---|
| 1350 | IF (debug) THEN |
---|
| 1351 | WRITE(numout,*) & |
---|
| 1352 | & ' OUT rest : z0, albedoVIS, albedoNIR, emis = ', & |
---|
| 1353 | & z0(itest,jtest),albedo(itest,jtest,1), & |
---|
| 1354 | & albedo(itest,jtest,2),emis(itest,jtest) |
---|
| 1355 | WRITE(numout,*) ' OUT rest : coastal and river flow = ', & |
---|
| 1356 | & coastalflow(itest,jtest), riverflow(itest,jtest) |
---|
| 1357 | WRITE(numout,*) ' OUT rest : tsol_rad, vevapp = ', & |
---|
| 1358 | & tsol_rad(itest,jtest), vevapp(itest,jtest) |
---|
| 1359 | WRITE(numout,*) ' OUT rest : temp_sol_new =', & |
---|
| 1360 | & temp_sol_NEW(itest,jtest) |
---|
| 1361 | ENDIF |
---|
| 1362 | |
---|
| 1363 | CALL barrier_para |
---|
| 1364 | CALL start_timer(timer_global) |
---|
| 1365 | CALL start_timer(timer_mpi) |
---|
| 1366 | |
---|
| 1367 | ENDIF |
---|
| 1368 | !----- |
---|
| 1369 | !---- Calling SECHIBA and doing the number crunching. |
---|
| 1370 | !---- Note that for the first time step SECHIBA is called twice. |
---|
| 1371 | !---- |
---|
| 1372 | !---- All H_2O fluxes are now in Kg/m^2s |
---|
| 1373 | !----- |
---|
| 1374 | IF (longprint) THEN |
---|
| 1375 | WRITE(numout,*) "dim2_driver ",it_force |
---|
| 1376 | WRITE(numout,*) ">> Index of land points =",kindex |
---|
| 1377 | WRITE(numout,*) "Lowest level wind speed North = ", & |
---|
| 1378 | & (/ ( for_u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1379 | WRITE(numout,*) "Lowest level wind speed East = ", & |
---|
| 1380 | & (/ ( for_v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1381 | WRITE(numout,*) "z0 ; Surface roughness = ", & |
---|
| 1382 | & (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1383 | WRITE(numout,*) "Height of first layer = ", & |
---|
| 1384 | & (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1385 | WRITE(numout,*) "Lowest level specific humidity = ", & |
---|
| 1386 | & (/ ( for_qair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1387 | WRITE(numout,*) "Rain precipitation = ", & |
---|
| 1388 | & (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /) |
---|
| 1389 | WRITE(numout,*) "Snow precipitation = ", & |
---|
| 1390 | & (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /) |
---|
| 1391 | WRITE(numout,*) "Down-welling long-wave flux = ", & |
---|
| 1392 | & (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1393 | WRITE(numout,*) "Net surface short-wave flux = ", & |
---|
| 1394 | & (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1395 | WRITE(numout,*) "Downwelling surface short-wave flux = ", & |
---|
| 1396 | & (/ ( for_swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1397 | WRITE(numout,*) "Air temperature in Kelvin = ", & |
---|
| 1398 | & (/ ( for_tair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1399 | WRITE(numout,*) "Air potential energy = ", & |
---|
| 1400 | & (/ ( for_eair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1401 | WRITE(numout,*) "CO2 concentration in the canopy = ", & |
---|
| 1402 | & (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1403 | WRITE(numout,*) "Coeficients A from the PBL resolution = ", & |
---|
| 1404 | & (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1405 | WRITE(numout,*) "One for T and another for q = ", & |
---|
| 1406 | & (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1407 | WRITE(numout,*) "Coeficients B from the PBL resolution = ", & |
---|
| 1408 | & (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1409 | WRITE(numout,*) "One for T and another for q = ", & |
---|
| 1410 | & (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1411 | WRITE(numout,*) "Cdrag = ", & |
---|
| 1412 | & (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1413 | WRITE(numout,*) "Lowest level pressure = ", & |
---|
| 1414 | & (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1415 | WRITE(numout,*) "Geographical coordinates lon = ", & |
---|
| 1416 | & (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1417 | WRITE(numout,*) "Geographical coordinates lat = ", & |
---|
| 1418 | & (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1419 | WRITE(numout,*) "Fraction of continent in the grid = ", & |
---|
| 1420 | & (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1421 | ENDIF |
---|
| 1422 | |
---|
| 1423 | CALL intersurf_main & |
---|
| 1424 | & (istp, iim, jjm, nbindex, kindex, dt, & |
---|
| 1425 | & first_CALL, last_CALL, lon, lat, for_contfrac, for_neighbours, for_resolution, date0, & |
---|
| 1426 | ! first level conditions |
---|
| 1427 | & zlev_vec, for_u, for_v, & |
---|
| 1428 | & for_qair, for_tair, for_eair, for_ccanopy, & |
---|
| 1429 | ! Variables for the implicit coupling |
---|
| 1430 | & cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, & |
---|
| 1431 | ! Rain, snow, radiation and surface pressure |
---|
| 1432 | & precip_rain, precip_snow, & |
---|
| 1433 | & for_lwdown, for_swnet, for_swdown, for_psurf, & |
---|
| 1434 | ! Output : Fluxes |
---|
| 1435 | & vevapp, fluxsens, fluxlat, coastalflow, riverflow, & |
---|
| 1436 | ! Surface temperatures and surface properties |
---|
| 1437 | & tsol_rad, temp_sol_NEW, qsurf, albedo, emis, z0) |
---|
| 1438 | !------- |
---|
| 1439 | IF (longprint) THEN |
---|
| 1440 | WRITE(numout,*) "dim2_driver outputs" |
---|
| 1441 | ! Output : Fluxes |
---|
| 1442 | WRITE(numout,*) "vevapp ; Total of evaporation = ", & |
---|
| 1443 | & (/ ( vevapp(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1444 | WRITE(numout,*) "Sensible heat flux = ", & |
---|
| 1445 | & (/ ( fluxsens(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1446 | WRITE(numout,*) "Latent heat flux = ", & |
---|
| 1447 | & (/ ( fluxlat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1448 | WRITE(numout,*) "coastalflow ; Diffuse flow of water into the ocean (m^3/dt) = ", & |
---|
| 1449 | & (/ ( coastalflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1450 | WRITE(numout,*) "riverflow ; Largest rivers flowing into the ocean (m^3/dt) = ", & |
---|
| 1451 | & (/ ( riverflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1452 | ! Surface temperatures and surface properties |
---|
| 1453 | WRITE(numout,*) "tsol_rad ; Radiative surface temperature = ", & |
---|
| 1454 | & (/ ( tsol_rad(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1455 | WRITE(numout,*) "temp_sol_new ; New soil temperature = ", & |
---|
| 1456 | & (/ ( temp_sol_NEW(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1457 | WRITE(numout,*) "qsurf ; Surface specific humidity = ", & |
---|
| 1458 | & (/ ( qsurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1459 | WRITE(numout,*) "albedoVIS = ", & |
---|
| 1460 | & (/ ( albedo(ilandindex(ik), jlandindex(ik), 1),ik=1,nbindex ) /) |
---|
| 1461 | WRITE(numout,*) "albedoNIR = ", & |
---|
| 1462 | & (/ ( albedo(ilandindex(ik), jlandindex(ik), 2),ik=1,nbindex ) /) |
---|
| 1463 | WRITE(numout,*) "emis ; Emissivity = ", & |
---|
| 1464 | & (/ ( emis(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1465 | WRITE(numout,*) "z0 ; Surface roughness = ", & |
---|
| 1466 | & (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) |
---|
| 1467 | ENDIF |
---|
| 1468 | !----- |
---|
| 1469 | dtdt(:,:) = zero |
---|
| 1470 | DO ik=1,nbindex |
---|
| 1471 | i=ilandindex(ik) |
---|
| 1472 | j=jlandindex(ik) |
---|
| 1473 | dtdt(i,j) = ABS(temp_sol_NEW(i,j)-temp_sol_old(i,j))/dt |
---|
| 1474 | ENDDO |
---|
| 1475 | !----- |
---|
| 1476 | !---- Test if the point with the largest change has more than 5K per dt |
---|
| 1477 | !----- |
---|
| 1478 | IF (MAXVAL(dtdt(:,:)) > 5./dt) THEN |
---|
| 1479 | ml = MAXLOC(dtdt) |
---|
| 1480 | CALL ju2ymds(julian, yy, mm, dd, ss) |
---|
| 1481 | WRITE(numout,"('ATT :',I5,' big temperature jumps on ', & |
---|
| 1482 | & I4,'-',I2.2,'-',I2.2,':',F8.4)") & |
---|
| 1483 | & COUNT(dtdt(:,:) > 5./dt),yy,mm,dd,ss/3600. |
---|
| 1484 | WRITE(numout,*) & |
---|
| 1485 | & 'Maximum change of surface temperature located at :', & |
---|
| 1486 | & lon(ml(1),ml(2)),lat(ml(1),ml(2)) |
---|
| 1487 | WRITE(numout,*) 'Coordinates in grid space: ',ml(1),ml(2) |
---|
| 1488 | WRITE(numout,*) 'Change from ',temp_sol_old(ml(1),ml(2)), & |
---|
| 1489 | & ' to ',temp_sol_new(ml(1),ml(2)),& |
---|
| 1490 | & 'with sw_in = ',for_swnet(ml(1),ml(2)) |
---|
| 1491 | old_tair = & |
---|
| 1492 | & (old_eair(ml(1),ml(2))-cte_grav*old_zlev(ml(1),ml(2)))/cp_air |
---|
| 1493 | WRITE(numout,*) 'Air temperature change from ',old_tair, & |
---|
| 1494 | & ' to ',for_tair(ml(1),ml(2)) |
---|
| 1495 | WRITE(numout,*) 'Max of dtdt : ',dtdt(ml(1),ml(2)),' with dt = ',dt |
---|
| 1496 | ENDIF |
---|
| 1497 | !----- |
---|
| 1498 | temp_sol_old(:,:) = temp_sol_NEW(:,:) |
---|
| 1499 | !----- |
---|
| 1500 | !---- PRINT output value for debug |
---|
| 1501 | !----- |
---|
| 1502 | IF (debug) THEN |
---|
| 1503 | WRITE(numout,*) ' OUT : z0, albedoVIS, albedoNIR, emis = ', & |
---|
| 1504 | & z0(itest,jtest),albedo(itest,jtest,1), & |
---|
| 1505 | & albedo(itest,jtest,2),emis(itest,jtest) |
---|
| 1506 | WRITE(numout,*) ' OUT : coastal and river flow = ',& |
---|
| 1507 | & coastalflow(itest,jtest), riverflow(itest,jtest) |
---|
| 1508 | WRITE(numout,*) ' OUT : tsol_rad, vevapp = ', & |
---|
| 1509 | & tsol_rad(itest,jtest), vevapp(itest,jtest) |
---|
| 1510 | WRITE(numout,*) ' OUT : temp_sol_new =', temp_sol_NEW(itest,jtest) |
---|
| 1511 | ENDIF |
---|
| 1512 | !----- |
---|
| 1513 | !---- Give some variables to the output package |
---|
| 1514 | !---- for saving on the history tape |
---|
| 1515 | !----- |
---|
| 1516 | IF (debug) WRITE(numout,*) 'history written for ', istp |
---|
| 1517 | !----- |
---|
| 1518 | istp_old = istp |
---|
| 1519 | istp = istp+1 |
---|
| 1520 | !----- |
---|
| 1521 | old_zlev(:,:) = zlev_vec(:,:) |
---|
| 1522 | old_qair(:,:) = for_qair(:,:) |
---|
| 1523 | old_eair(:,:) = for_eair(:,:) |
---|
| 1524 | !----- |
---|
| 1525 | is = is + 1 |
---|
| 1526 | ENDDO |
---|
| 1527 | split_start = 1 |
---|
| 1528 | !!$ it_force =it_force+1 |
---|
| 1529 | IF (it==itau_fin-1) THEN |
---|
| 1530 | CALL stop_timer(timer_global) |
---|
| 1531 | CALL stop_timer(timer_mpi) |
---|
| 1532 | ENDIF |
---|
| 1533 | it = it + 1 |
---|
| 1534 | !!! CALL histsync |
---|
| 1535 | ENDDO |
---|
| 1536 | !- |
---|
| 1537 | ! Archive in restart file the prognostic variables |
---|
| 1538 | !- |
---|
| 1539 | IF (debug) WRITE(numout,*) 'Write the restart for the driver', istp_old |
---|
| 1540 | !- |
---|
| 1541 | var_name = 'fluxsens' |
---|
| 1542 | IF (is_root_prc) & |
---|
| 1543 | ALLOCATE(fluxsens_g(iim_g,jjm_g)) |
---|
| 1544 | CALL gather2D(fluxsens , fluxsens_g) |
---|
| 1545 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, fluxsens_g) |
---|
| 1546 | IF (is_root_prc) & |
---|
| 1547 | DEALLOCATE(fluxsens_g) |
---|
| 1548 | |
---|
| 1549 | var_name = 'vevapp' |
---|
| 1550 | IF (is_root_prc) & |
---|
| 1551 | ALLOCATE(vevapp_g(iim_g,jjm_g)) |
---|
| 1552 | CALL gather2D( vevapp, vevapp_g) |
---|
| 1553 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, vevapp_g) |
---|
| 1554 | IF (is_root_prc) & |
---|
| 1555 | DEALLOCATE(vevapp_g) |
---|
| 1556 | |
---|
| 1557 | var_name = 'zlev_old' |
---|
| 1558 | IF (is_root_prc) & |
---|
| 1559 | ALLOCATE(old_zlev_g(iim_g,jjm_g)) |
---|
| 1560 | CALL gather2D( old_zlev, old_zlev_g) |
---|
| 1561 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_zlev_g) |
---|
| 1562 | IF (is_root_prc) & |
---|
| 1563 | DEALLOCATE(old_zlev_g) |
---|
| 1564 | |
---|
| 1565 | var_name = 'qair_old' |
---|
| 1566 | IF (is_root_prc) & |
---|
| 1567 | ALLOCATE(old_qair_g(iim_g,jjm_g)) |
---|
| 1568 | CALL gather2D( old_qair, old_qair_g) |
---|
| 1569 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_qair_g) |
---|
| 1570 | IF (is_root_prc) & |
---|
| 1571 | DEALLOCATE(old_qair_g) |
---|
| 1572 | |
---|
| 1573 | var_name = 'eair_old' |
---|
| 1574 | IF (is_root_prc) & |
---|
| 1575 | ALLOCATE(old_eair_g(iim_g,jjm_g)) |
---|
| 1576 | CALL gather2D( old_eair, old_eair_g) |
---|
| 1577 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_eair_g) |
---|
| 1578 | IF (is_root_prc) & |
---|
| 1579 | DEALLOCATE(old_eair_g) |
---|
| 1580 | |
---|
| 1581 | var_name = 'rau_old' |
---|
| 1582 | IF (is_root_prc) & |
---|
| 1583 | ALLOCATE(for_rau_g(iim_g,jjm_g)) |
---|
| 1584 | CALL gather2D( for_rau, for_rau_g) |
---|
| 1585 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, for_rau_g) |
---|
| 1586 | IF (is_root_prc) & |
---|
| 1587 | DEALLOCATE(for_rau_g) |
---|
| 1588 | |
---|
| 1589 | IF (is_root_prc) & |
---|
| 1590 | ALLOCATE(albedo_g(iim_g,jjm_g)) |
---|
| 1591 | var_name= 'albedo_vis' |
---|
| 1592 | albedo_vis(:,:)=albedo(:,:,1) |
---|
| 1593 | CALL gather2D(albedo_vis,albedo_g) |
---|
| 1594 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, albedo_g) |
---|
| 1595 | ! |
---|
| 1596 | var_name= 'albedo_nir' |
---|
| 1597 | albedo_nir(:,:)=albedo(:,:,2) |
---|
| 1598 | CALL gather2D(albedo_nir,albedo_g) |
---|
| 1599 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, albedo_g) |
---|
| 1600 | IF (is_root_prc) & |
---|
| 1601 | DEALLOCATE(albedo_g) |
---|
| 1602 | |
---|
| 1603 | IF (is_root_prc) & |
---|
| 1604 | ALLOCATE(z0_g(iim_g,jjm_g)) |
---|
| 1605 | var_name= 'z0' |
---|
| 1606 | CALL gather2D(z0,z0_g) |
---|
| 1607 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, z0_g) |
---|
| 1608 | IF (is_root_prc) & |
---|
| 1609 | DEALLOCATE(z0_g) |
---|
| 1610 | |
---|
| 1611 | if (.NOT. is_watchout) THEN |
---|
| 1612 | var_name = 'petAcoef' |
---|
| 1613 | IF (is_root_prc) & |
---|
| 1614 | ALLOCATE(petAcoef_g(iim_g,jjm_g)) |
---|
| 1615 | CALL gather2D( petAcoef, petAcoef_g) |
---|
| 1616 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, petAcoef_g) |
---|
| 1617 | IF (is_root_prc) & |
---|
| 1618 | DEALLOCATE(petAcoef_g) |
---|
| 1619 | |
---|
| 1620 | var_name = 'petBcoef' |
---|
| 1621 | IF (is_root_prc) & |
---|
| 1622 | ALLOCATE(petBcoef_g(iim_g,jjm_g)) |
---|
| 1623 | CALL gather2D( petBcoef, petBcoef_g) |
---|
| 1624 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, petBcoef_g) |
---|
| 1625 | IF (is_root_prc) & |
---|
| 1626 | DEALLOCATE(petBcoef_g) |
---|
| 1627 | |
---|
| 1628 | var_name = 'peqAcoef' |
---|
| 1629 | IF (is_root_prc) & |
---|
| 1630 | ALLOCATE(peqAcoef_g(iim_g,jjm_g)) |
---|
| 1631 | CALL gather2D( peqAcoef, peqAcoef_g) |
---|
| 1632 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, peqAcoef_g) |
---|
| 1633 | IF (is_root_prc) & |
---|
| 1634 | DEALLOCATE(peqAcoef_g) |
---|
| 1635 | |
---|
| 1636 | var_name = 'peqBcoef' |
---|
| 1637 | IF (is_root_prc) & |
---|
| 1638 | ALLOCATE(peqBcoef_g(iim_g,jjm_g)) |
---|
| 1639 | CALL gather2D( peqBcoef, peqBcoef_g) |
---|
| 1640 | IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, peqBcoef_g) |
---|
| 1641 | IF (is_root_prc) & |
---|
| 1642 | DEALLOCATE(peqBcoef_g) |
---|
| 1643 | ENDIF |
---|
| 1644 | !- |
---|
| 1645 | IF (debug) WRITE(numout,*) 'Restart for the driver written' |
---|
| 1646 | !===================================================================== |
---|
| 1647 | !- 5.0 Closing all files |
---|
| 1648 | !===================================================================== |
---|
| 1649 | CALL flinclo(force_id) |
---|
| 1650 | IF ( debug ) WRITE(numout,*) 'FLIN CLOSED' |
---|
| 1651 | CALL histclo |
---|
| 1652 | IF ( debug ) WRITE(numout,*) 'HIST CLOSED' |
---|
| 1653 | |
---|
| 1654 | IF(is_root_prc) THEN |
---|
| 1655 | CALL restclo |
---|
| 1656 | IF ( debug ) WRITE(numout,*) 'REST CLOSED' |
---|
| 1657 | CALL getin_dump |
---|
| 1658 | IF ( debug ) WRITE(numout,*) 'GETIN CLOSED' |
---|
| 1659 | ENDIF |
---|
| 1660 | !--------------- |
---|
| 1661 | CALL finalize_para(timer_global,timer_mpi) |
---|
| 1662 | !- |
---|
| 1663 | STOP 'END of dim2_driver' |
---|
| 1664 | !--------------- |
---|
| 1665 | END PROGRAM driver |
---|