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