[222] | 1 | C====================================================================== |
---|
| 2 | PROGRAM newstart |
---|
| 3 | c======================================================================= |
---|
| 4 | c |
---|
| 5 | c |
---|
| 6 | c Auteur: Christophe Hourdin/Francois Forget/Yann Wanherdrick |
---|
| 7 | c ------ |
---|
| 8 | c Derniere modif : 12/03 |
---|
| 9 | c |
---|
| 10 | c |
---|
| 11 | c Objet: Create or modify the initial state for the LMD Mars GCM |
---|
| 12 | c ----- (fichiers NetCDF start et startfi) |
---|
| 13 | c |
---|
| 14 | c |
---|
| 15 | c======================================================================= |
---|
| 16 | |
---|
| 17 | use infotrac, only: iniadvtrac, nqtot, tname |
---|
| 18 | USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice |
---|
| 19 | USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat |
---|
| 20 | USE surfdat_h, ONLY: phisfi, albedodat, |
---|
| 21 | & zmea, zstd, zsig, zgam, zthe |
---|
| 22 | USE comgeomfi_h, ONLY: lati, long, area |
---|
| 23 | use datafile_mod, only: datadir |
---|
| 24 | ! to use 'getin' |
---|
[227] | 25 | ! USE ioipsl_getincom, only: getin |
---|
| 26 | USE ioipsl_getincom_p, only: getin_p |
---|
[222] | 27 | use control_mod, only: day_step, iphysiq, anneeref |
---|
| 28 | use phyredem, only: physdem0, physdem1 |
---|
| 29 | use iostart, only: open_startphy |
---|
| 30 | use comgeomphy, only: initcomgeomphy |
---|
| 31 | use slab_ice_h, only:noceanmx |
---|
| 32 | implicit none |
---|
| 33 | |
---|
| 34 | #include "dimensions.h" |
---|
[227] | 35 | !#include "dimphys.h" |
---|
| 36 | integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm) |
---|
| 37 | !#include "planete.h" |
---|
[222] | 38 | #include "paramet.h" |
---|
| 39 | #include "comconst.h" |
---|
| 40 | #include "comvert.h" |
---|
| 41 | #include "comgeom2.h" |
---|
| 42 | !#include "control.h" |
---|
| 43 | #include "logic.h" |
---|
| 44 | #include "ener.h" |
---|
| 45 | #include "temps.h" |
---|
| 46 | #include "comdissnew.h" |
---|
| 47 | #include "serre.h" |
---|
| 48 | #include "netcdf.inc" |
---|
| 49 | !#include "advtrac.h" |
---|
| 50 | c======================================================================= |
---|
| 51 | c Declarations |
---|
| 52 | c======================================================================= |
---|
| 53 | |
---|
| 54 | c Variables dimension du fichier "start_archive" |
---|
| 55 | c------------------------------------ |
---|
| 56 | CHARACTER relief*3 |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | c Variables pour les lectures NetCDF des fichiers "start_archive" |
---|
| 60 | c-------------------------------------------------- |
---|
| 61 | INTEGER nid_dyn, nid_fi,nid,nvarid |
---|
| 62 | INTEGER length |
---|
| 63 | parameter (length = 100) |
---|
| 64 | INTEGER tab0 |
---|
| 65 | INTEGER NB_ETATMAX |
---|
| 66 | parameter (NB_ETATMAX = 100) |
---|
| 67 | |
---|
| 68 | REAL date |
---|
| 69 | REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec |
---|
| 70 | |
---|
| 71 | c Variable histoire |
---|
| 72 | c------------------ |
---|
| 73 | REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants |
---|
| 74 | REAL phis(iip1,jjp1) |
---|
| 75 | REAL,ALLOCATABLE :: q(:,:,:,:) ! champs advectes |
---|
| 76 | |
---|
| 77 | c autre variables dynamique nouvelle grille |
---|
| 78 | c------------------------------------------ |
---|
| 79 | REAL pks(iip1,jjp1) |
---|
| 80 | REAL w(iip1,jjp1,llm+1) |
---|
| 81 | REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) |
---|
| 82 | ! REAL dv(ip1jm,llm),du(ip1jmp1,llm) |
---|
| 83 | ! REAL dh(ip1jmp1,llm),dp(ip1jmp1) |
---|
| 84 | REAL phi(iip1,jjp1,llm) |
---|
| 85 | |
---|
| 86 | integer klatdat,klongdat |
---|
| 87 | PARAMETER (klatdat=180,klongdat=360) |
---|
| 88 | |
---|
| 89 | c Physique sur grille scalaire |
---|
| 90 | c---------------------------- |
---|
| 91 | real zmeaS(iip1,jjp1),zstdS(iip1,jjp1) |
---|
| 92 | real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1) |
---|
| 93 | |
---|
| 94 | c variable physique |
---|
| 95 | c------------------ |
---|
| 96 | REAL tsurf(ngridmx) ! surface temperature |
---|
| 97 | REAL tsoil(ngridmx,nsoilmx) ! soil temperature |
---|
| 98 | ! REAL co2ice(ngridmx) ! CO2 ice layer |
---|
| 99 | REAL emis(ngridmx) ! surface emissivity |
---|
| 100 | real emisread ! added by RW |
---|
| 101 | REAL,ALLOCATABLE :: qsurf(:,:) |
---|
[227] | 102 | REAL q2(ngridmx,llm+1) |
---|
[222] | 103 | ! REAL rnaturfi(ngridmx) |
---|
| 104 | real alb(iip1,jjp1),albfi(ngridmx) ! albedos |
---|
| 105 | real ith(iip1,jjp1,nsoilmx),ithfi(ngridmx,nsoilmx) ! thermal inertia (3D) |
---|
| 106 | real surfith(iip1,jjp1),surfithfi(ngridmx) ! surface thermal inertia (2D) |
---|
| 107 | REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx) |
---|
| 108 | |
---|
| 109 | INTEGER i,j,l,isoil,ig,idum |
---|
| 110 | real mugaz ! molar mass of the atmosphere |
---|
| 111 | |
---|
| 112 | integer ierr |
---|
| 113 | |
---|
| 114 | REAL rnat(ngridmx) |
---|
| 115 | REAL tslab(ngridmx,nsoilmx) ! slab ocean temperature |
---|
| 116 | REAL pctsrf_sic(ngridmx) ! sea ice cover |
---|
| 117 | REAL tsea_ice(ngridmx) ! temperature sea_ice |
---|
| 118 | REAL sea_ice(ngridmx) ! mass of sea ice (kg/m2) |
---|
| 119 | |
---|
| 120 | c Variables on the new grid along scalar points |
---|
| 121 | c------------------------------------------------------ |
---|
| 122 | ! REAL p(iip1,jjp1) |
---|
| 123 | REAL t(iip1,jjp1,llm) |
---|
| 124 | REAL tset(iip1,jjp1,llm) |
---|
| 125 | real phisold_newgrid(iip1,jjp1) |
---|
| 126 | REAL :: teta(iip1, jjp1, llm) |
---|
| 127 | REAL :: pk(iip1,jjp1,llm) |
---|
| 128 | REAL :: pkf(iip1,jjp1,llm) |
---|
| 129 | REAL :: ps(iip1, jjp1) |
---|
| 130 | REAL :: masse(iip1,jjp1,llm) |
---|
| 131 | REAL :: xpn,xps,xppn(iim),xpps(iim) |
---|
| 132 | REAL :: p3d(iip1, jjp1, llm+1) |
---|
| 133 | REAL :: beta(iip1,jjp1,llm) |
---|
| 134 | ! REAL dteta(ip1jmp1,llm) |
---|
| 135 | |
---|
| 136 | c Variable de l'ancienne grille |
---|
| 137 | c------------------------------ |
---|
| 138 | real time |
---|
| 139 | real tab_cntrl(100) |
---|
| 140 | real tab_cntrl_bis(100) |
---|
| 141 | |
---|
| 142 | c variables diverses |
---|
| 143 | c------------------- |
---|
| 144 | real choix_1,pp |
---|
| 145 | character*80 fichnom |
---|
| 146 | character*250 filestring |
---|
| 147 | integer Lmodif,iq,thermo |
---|
| 148 | character modif*20 |
---|
| 149 | real z_reel(iip1,jjp1) |
---|
| 150 | real tsud,albsud,alb_bb,ith_bb,Tiso,Tabove |
---|
| 151 | real ptoto,pcap,patm,airetot,ptotn,patmn,psea |
---|
| 152 | ! real ssum |
---|
| 153 | character*1 yes |
---|
| 154 | logical :: flagtset=.false. , flagps0=.false. |
---|
| 155 | real val, val2, val3 ! to store temporary variables |
---|
| 156 | real :: iceith=2000 ! thermal inertia of subterranean ice |
---|
| 157 | integer iref,jref |
---|
| 158 | |
---|
| 159 | INTEGER :: itau |
---|
| 160 | |
---|
| 161 | INTEGER :: nq,numvanle |
---|
| 162 | character(len=20) :: txt ! to store some text |
---|
| 163 | character(len=50) :: surfacefile ! "surface.nc" (or equivalent file) |
---|
| 164 | character(len=150) :: longline |
---|
| 165 | integer :: count |
---|
| 166 | real :: profile(llm+1) ! to store an atmospheric profile + surface value |
---|
| 167 | |
---|
| 168 | ! added by RW for test |
---|
| 169 | real pmean, phi0 |
---|
| 170 | |
---|
| 171 | ! added by BC for equilibrium temperature startup |
---|
| 172 | real teque |
---|
| 173 | |
---|
| 174 | ! added by BC for cloud fraction setup |
---|
[227] | 175 | REAL hice(ngridmx),cloudfrac(ngridmx,llm) |
---|
[222] | 176 | REAL totalfrac(ngridmx) |
---|
| 177 | |
---|
| 178 | |
---|
| 179 | ! added by RW for nuketharsis |
---|
| 180 | real fact1 |
---|
| 181 | real fact2 |
---|
| 182 | |
---|
| 183 | |
---|
| 184 | c sortie visu pour les champs dynamiques |
---|
| 185 | c--------------------------------------- |
---|
| 186 | ! INTEGER :: visuid |
---|
| 187 | ! real :: time_step,t_ops,t_wrt |
---|
| 188 | ! CHARACTER*80 :: visu_file |
---|
| 189 | |
---|
| 190 | cpp = 0. |
---|
| 191 | preff = 0. |
---|
| 192 | pa = 0. ! to ensure disaster rather than minor error if we don`t |
---|
| 193 | ! make deliberate choice of these values elsewhere. |
---|
| 194 | |
---|
| 195 | ! initialize "serial/parallel" related stuff |
---|
| 196 | CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) |
---|
| 197 | call initcomgeomphy |
---|
| 198 | ! Load tracer number and names: |
---|
| 199 | call iniadvtrac(nqtot,numvanle) |
---|
| 200 | ! allocate arrays |
---|
| 201 | allocate(q(iip1,jjp1,llm,nqtot)) |
---|
| 202 | allocate(qsurf(ngridmx,nqtot)) |
---|
| 203 | |
---|
| 204 | c======================================================================= |
---|
| 205 | c Choice of the start file(s) to use |
---|
| 206 | c======================================================================= |
---|
| 207 | write(*,*) 'From which kind of files do you want to create new', |
---|
| 208 | . 'start and startfi files' |
---|
| 209 | write(*,*) ' 0 - from a file start_archive' |
---|
| 210 | write(*,*) ' 1 - from files start and startfi' |
---|
| 211 | |
---|
| 212 | c----------------------------------------------------------------------- |
---|
| 213 | c Open file(s) to modify (start or start_archive) |
---|
| 214 | c----------------------------------------------------------------------- |
---|
| 215 | |
---|
| 216 | DO |
---|
| 217 | read(*,*,iostat=ierr) choix_1 |
---|
| 218 | if ((choix_1 /= 0).OR.(choix_1 /=1)) EXIT |
---|
| 219 | ENDDO |
---|
| 220 | |
---|
| 221 | c Open start_archive |
---|
| 222 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 223 | if (choix_1.eq.0) then |
---|
| 224 | |
---|
| 225 | write(*,*) 'Creating start files from:' |
---|
| 226 | write(*,*) './start_archive.nc' |
---|
| 227 | write(*,*) |
---|
| 228 | fichnom = 'start_archive.nc' |
---|
| 229 | ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) |
---|
| 230 | IF (ierr.NE.NF_NOERR) THEN |
---|
| 231 | write(6,*)' Problem opening file:',fichnom |
---|
| 232 | write(6,*)' ierr = ', ierr |
---|
| 233 | CALL ABORT |
---|
| 234 | ENDIF |
---|
| 235 | tab0 = 50 |
---|
| 236 | Lmodif = 1 |
---|
| 237 | |
---|
| 238 | c OR open start and startfi files |
---|
| 239 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 240 | else |
---|
| 241 | write(*,*) 'Creating start files from:' |
---|
| 242 | write(*,*) './start.nc and ./startfi.nc' |
---|
| 243 | write(*,*) |
---|
| 244 | fichnom = 'start.nc' |
---|
| 245 | ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_dyn) |
---|
| 246 | IF (ierr.NE.NF_NOERR) THEN |
---|
| 247 | write(6,*)' Problem opening file:',fichnom |
---|
| 248 | write(6,*)' ierr = ', ierr |
---|
| 249 | CALL ABORT |
---|
| 250 | ENDIF |
---|
| 251 | |
---|
| 252 | fichnom = 'startfi.nc' |
---|
| 253 | ierr = NF_OPEN (fichnom, NF_NOWRITE,nid_fi) |
---|
| 254 | IF (ierr.NE.NF_NOERR) THEN |
---|
| 255 | write(6,*)' Problem opening file:',fichnom |
---|
| 256 | write(6,*)' ierr = ', ierr |
---|
| 257 | CALL ABORT |
---|
| 258 | ENDIF |
---|
| 259 | |
---|
| 260 | tab0 = 0 |
---|
| 261 | Lmodif = 0 |
---|
| 262 | |
---|
| 263 | endif |
---|
| 264 | |
---|
| 265 | |
---|
| 266 | c======================================================================= |
---|
| 267 | c INITIALISATIONS DIVERSES |
---|
| 268 | c======================================================================= |
---|
| 269 | |
---|
| 270 | ! Initialize global tracer indexes (stored in tracer.h) |
---|
| 271 | ! ... this has to be done before phyetat0 |
---|
| 272 | call initracer(ngridmx,nqtot,tname) |
---|
| 273 | |
---|
| 274 | ! Take care of arrays in common modules |
---|
| 275 | ! ALLOCATE ARRAYS in surfdat_h (if not already done, e.g. when using start_archive) |
---|
| 276 | IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngridmx)) |
---|
| 277 | IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngridmx)) |
---|
| 278 | IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngridmx)) |
---|
| 279 | IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngridmx)) |
---|
| 280 | IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngridmx)) |
---|
| 281 | IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngridmx)) |
---|
| 282 | IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngridmx)) |
---|
| 283 | ! ALLOCATE ARRAYS in comsoil_h (if not already done) |
---|
| 284 | IF (.not.ALLOCATED(layer)) |
---|
| 285 | . ALLOCATE(layer(nsoilmx)) |
---|
| 286 | IF (.not.ALLOCATED(mlayer)) |
---|
| 287 | . ALLOCATE(mlayer(0:nsoilmx-1)) |
---|
| 288 | IF (.not.ALLOCATED(inertiedat)) |
---|
| 289 | . ALLOCATE(inertiedat(ngridmx,nsoilmx)) |
---|
| 290 | ! ALLOCATE ARRAYS IN comgeomfi_h (done in inifis usually) |
---|
| 291 | IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngridmx)) |
---|
| 292 | IF (.not. ALLOCATED(long)) ALLOCATE(long(ngridmx)) |
---|
| 293 | IF (.not. ALLOCATED(area)) ALLOCATE(area(ngridmx)) |
---|
| 294 | |
---|
| 295 | c----------------------------------------------------------------------- |
---|
| 296 | c Lecture du tableau des parametres du run (pour la dynamique) |
---|
| 297 | c----------------------------------------------------------------------- |
---|
| 298 | if (choix_1.eq.0) then |
---|
| 299 | |
---|
| 300 | write(*,*) 'reading tab_cntrl START_ARCHIVE' |
---|
| 301 | c |
---|
| 302 | ierr = NF_INQ_VARID (nid, "controle", nvarid) |
---|
| 303 | #ifdef NC_DOUBLE |
---|
| 304 | ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) |
---|
| 305 | #else |
---|
| 306 | ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) |
---|
| 307 | #endif |
---|
| 308 | c |
---|
| 309 | else if (choix_1.eq.1) then |
---|
| 310 | |
---|
| 311 | write(*,*) 'reading tab_cntrl START' |
---|
| 312 | c |
---|
| 313 | ierr = NF_INQ_VARID (nid_dyn, "controle", nvarid) |
---|
| 314 | #ifdef NC_DOUBLE |
---|
| 315 | ierr = NF_GET_VAR_DOUBLE(nid_dyn, nvarid, tab_cntrl) |
---|
| 316 | #else |
---|
| 317 | ierr = NF_GET_VAR_REAL(nid_dyn, nvarid, tab_cntrl) |
---|
| 318 | #endif |
---|
| 319 | c |
---|
| 320 | write(*,*) 'reading tab_cntrl STARTFI' |
---|
| 321 | c |
---|
| 322 | ierr = NF_INQ_VARID (nid_fi, "controle", nvarid) |
---|
| 323 | #ifdef NC_DOUBLE |
---|
| 324 | ierr = NF_GET_VAR_DOUBLE(nid_fi, nvarid, tab_cntrl_bis) |
---|
| 325 | #else |
---|
| 326 | ierr = NF_GET_VAR_REAL(nid_fi, nvarid, tab_cntrl_bis) |
---|
| 327 | #endif |
---|
| 328 | c |
---|
| 329 | do i=1,50 |
---|
| 330 | tab_cntrl(i+50)=tab_cntrl_bis(i) |
---|
| 331 | enddo |
---|
| 332 | write(*,*) 'printing tab_cntrl', tab_cntrl |
---|
| 333 | do i=1,100 |
---|
| 334 | write(*,*) i,tab_cntrl(i) |
---|
| 335 | enddo |
---|
| 336 | |
---|
| 337 | ! Lmodif set to 0 to disable modifications possibility in phyeta0 |
---|
| 338 | write(*,*) 'Reading file START' |
---|
| 339 | fichnom = 'start.nc' |
---|
| 340 | CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse, |
---|
| 341 | . ps,phis,time) |
---|
| 342 | |
---|
| 343 | write(*,*) 'Reading file STARTFI' |
---|
| 344 | fichnom = 'startfi.nc' |
---|
[227] | 345 | CALL phyetat0 (ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx, |
---|
| 346 | . nqtot,day_ini,time, |
---|
[222] | 347 | . tsurf,tsoil,emis,q2,qsurf, !) ! temporary modif by RDW |
---|
| 348 | . cloudfrac,totalfrac,hice,rnat,pctsrf_sic,tslab,tsea_ice, |
---|
| 349 | . sea_ice) |
---|
| 350 | |
---|
| 351 | ! copy albedo and soil thermal inertia on (local) physics grid |
---|
| 352 | do i=1,ngridmx |
---|
| 353 | albfi(i) = albedodat(i) |
---|
| 354 | do j=1,nsoilmx |
---|
| 355 | ithfi(i,j) = inertiedat(i,j) |
---|
| 356 | enddo |
---|
| 357 | ! build a surfithfi(:) using 1st layer of ithfi(:), which might |
---|
| 358 | ! be needed later on if reinitializing soil thermal inertia |
---|
| 359 | surfithfi(i)=ithfi(i,1) |
---|
| 360 | enddo |
---|
| 361 | ! also copy albedo and soil thermal inertia on (local) dynamics grid |
---|
| 362 | ! so that options below can manipulate either (but must then ensure |
---|
| 363 | ! to correctly recast things on physics grid) |
---|
| 364 | call gr_fi_dyn(1,ngridmx,iip1,jjp1,albfi,alb) |
---|
| 365 | call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) |
---|
| 366 | call gr_fi_dyn(1,ngridmx,iip1,jjp1,surfithfi,surfith) |
---|
| 367 | |
---|
| 368 | endif |
---|
| 369 | c----------------------------------------------------------------------- |
---|
| 370 | c Initialisation des constantes dynamique |
---|
| 371 | c----------------------------------------------------------------------- |
---|
| 372 | |
---|
| 373 | kappa = tab_cntrl(9) |
---|
| 374 | etot0 = tab_cntrl(12) |
---|
| 375 | ptot0 = tab_cntrl(13) |
---|
| 376 | ztot0 = tab_cntrl(14) |
---|
| 377 | stot0 = tab_cntrl(15) |
---|
| 378 | ang0 = tab_cntrl(16) |
---|
| 379 | write(*,*) "Newstart: kappa,etot0,ptot0,ztot0,stot0,ang0" |
---|
| 380 | write(*,*) kappa,etot0,ptot0,ztot0,stot0,ang0 |
---|
| 381 | |
---|
| 382 | ! for vertical coordinate |
---|
| 383 | preff=tab_cntrl(18) ! reference surface pressure |
---|
| 384 | pa=tab_cntrl(17) ! reference pressure at which coord is purely pressure |
---|
| 385 | !NB: in start_archive files tab_cntrl(17)=tab_cntrl(18)=0 |
---|
| 386 | write(*,*) "Newstart: preff=",preff," pa=",pa |
---|
| 387 | yes=' ' |
---|
| 388 | do while ((yes.ne.'y').and.(yes.ne.'n')) |
---|
| 389 | write(*,*) "Change the values of preff and pa? (y/n)" |
---|
| 390 | read(*,fmt='(a)') yes |
---|
| 391 | end do |
---|
| 392 | if (yes.eq.'y') then |
---|
| 393 | write(*,*)"New value of reference surface pressure preff?" |
---|
| 394 | write(*,*)" (for Mars, typically preff=610)" |
---|
| 395 | read(*,*) preff |
---|
| 396 | write(*,*)"New value of reference pressure pa for purely" |
---|
| 397 | write(*,*)"pressure levels (for hybrid coordinates)?" |
---|
| 398 | write(*,*)" (for Mars, typically pa=20)" |
---|
| 399 | read(*,*) pa |
---|
| 400 | endif |
---|
| 401 | c----------------------------------------------------------------------- |
---|
| 402 | c Lecture du tab_cntrl et initialisation des constantes physiques |
---|
| 403 | c - pour start: Lmodif = 0 => pas de modifications possibles |
---|
| 404 | c (modif dans le tabfi de readfi + loin) |
---|
| 405 | c - pour start_archive: Lmodif = 1 => modifications possibles |
---|
| 406 | c----------------------------------------------------------------------- |
---|
| 407 | if (choix_1.eq.0) then |
---|
| 408 | ! tabfi requires that input file be first opened by open_startphy(fichnom) |
---|
| 409 | fichnom = 'start_archive.nc' |
---|
| 410 | call open_startphy(fichnom) |
---|
| 411 | call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad, |
---|
| 412 | . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) |
---|
| 413 | else if (choix_1.eq.1) then |
---|
| 414 | fichnom = 'startfi.nc' |
---|
| 415 | call open_startphy(fichnom) |
---|
| 416 | Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0 |
---|
| 417 | call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad, |
---|
| 418 | . p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time) |
---|
| 419 | endif |
---|
| 420 | |
---|
| 421 | rad = p_rad |
---|
| 422 | omeg = p_omeg |
---|
| 423 | g = p_g |
---|
| 424 | cpp = p_cpp |
---|
| 425 | mugaz = p_mugaz |
---|
| 426 | daysec = p_daysec |
---|
| 427 | |
---|
| 428 | kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW |
---|
| 429 | |
---|
| 430 | c======================================================================= |
---|
| 431 | c INITIALISATIONS DIVERSES |
---|
| 432 | c======================================================================= |
---|
| 433 | ! Load parameters from run.def file |
---|
| 434 | CALL defrun_new( 99, .TRUE. ) |
---|
| 435 | CALL iniconst |
---|
| 436 | CALL inigeom |
---|
| 437 | idum=-1 |
---|
| 438 | idum=0 |
---|
| 439 | |
---|
| 440 | c Initialisation coordonnees /aires |
---|
| 441 | c ------------------------------- |
---|
| 442 | ! Note: rlatu(:) and rlonv(:) are commons defined in "comgeom.h" |
---|
| 443 | ! rlatu() and rlonv() are given in radians |
---|
| 444 | latfi(1)=rlatu(1) |
---|
| 445 | lonfi(1)=0. |
---|
| 446 | DO j=2,jjm |
---|
| 447 | DO i=1,iim |
---|
| 448 | latfi((j-2)*iim+1+i)=rlatu(j) |
---|
| 449 | lonfi((j-2)*iim+1+i)=rlonv(i) |
---|
| 450 | ENDDO |
---|
| 451 | ENDDO |
---|
| 452 | latfi(ngridmx)=rlatu(jjp1) |
---|
| 453 | lonfi(ngridmx)=0. |
---|
| 454 | |
---|
| 455 | ! build airefi(), mesh area on physics grid |
---|
| 456 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) |
---|
| 457 | ! Poles are single points on physics grid |
---|
| 458 | airefi(1)=airefi(1)*iim |
---|
| 459 | airefi(ngridmx)=airefi(ngridmx)*iim |
---|
| 460 | |
---|
| 461 | ! also initialize various physics flags/settings which might be needed |
---|
| 462 | ! (for instance initracer needs to know about some flags, and/or |
---|
| 463 | ! 'datafile' path may be changed by user) |
---|
| 464 | call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys, |
---|
| 465 | & latfi,lonfi,airefi,rad,g,r,cpp) |
---|
| 466 | |
---|
| 467 | c======================================================================= |
---|
| 468 | c lecture topographie, albedo, inertie thermique, relief sous-maille |
---|
| 469 | c======================================================================= |
---|
| 470 | |
---|
| 471 | if (choix_1.eq.0) then ! for start_archive files, |
---|
| 472 | ! where an external "surface.nc" file is needed |
---|
| 473 | |
---|
| 474 | c do while((relief(1:3).ne.'mol').AND.(relief(1:3).ne.'pla')) |
---|
| 475 | c write(*,*) |
---|
| 476 | c write(*,*) 'choix du relief (mola,pla)' |
---|
| 477 | c write(*,*) '(Topographie MGS MOLA, plat)' |
---|
| 478 | c read(*,fmt='(a3)') relief |
---|
| 479 | relief="mola" |
---|
| 480 | c enddo |
---|
| 481 | |
---|
| 482 | ! First get the correct value of datadir, if not already done: |
---|
| 483 | ! default 'datadir' is set in "datafile_mod" |
---|
[227] | 484 | call getin_p("datadir",datadir) |
---|
[222] | 485 | write(*,*) 'Available surface data files are:' |
---|
| 486 | filestring='ls '//trim(datadir)//' | grep .nc' |
---|
| 487 | call system(filestring) |
---|
| 488 | |
---|
| 489 | write(*,*) '' |
---|
| 490 | write(*,*) 'Please choose the relevant file', |
---|
| 491 | & ' (e.g. "surface_mars.nc")' |
---|
| 492 | write(*,*) ' or "none" to not use any (and set planetary' |
---|
| 493 | write(*,*) ' albedo and surface thermal inertia)' |
---|
| 494 | read(*,fmt='(a50)') surfacefile |
---|
| 495 | |
---|
| 496 | if (surfacefile.ne."none") then |
---|
| 497 | |
---|
| 498 | CALL datareadnc(relief,surfacefile,phis,alb,surfith, |
---|
| 499 | & zmeaS,zstdS,zsigS,zgamS,ztheS) |
---|
| 500 | else |
---|
| 501 | ! specific case when not using a "surface.nc" file |
---|
| 502 | phis(:,:)=0 |
---|
| 503 | zmeaS(:,:)=0 |
---|
| 504 | zstdS(:,:)=0 |
---|
| 505 | zsigS(:,:)=0 |
---|
| 506 | zgamS(:,:)=0 |
---|
| 507 | ztheS(:,:)=0 |
---|
| 508 | |
---|
| 509 | write(*,*) "Enter value of albedo of the bare ground:" |
---|
| 510 | read(*,*) alb(1,1) |
---|
| 511 | alb(:,:)=alb(1,1) |
---|
| 512 | |
---|
| 513 | write(*,*) "Enter value of thermal inertia of soil:" |
---|
| 514 | read(*,*) surfith(1,1) |
---|
| 515 | surfith(:,:)=surfith(1,1) |
---|
| 516 | |
---|
| 517 | endif ! of if (surfacefile.ne."none") |
---|
| 518 | |
---|
| 519 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) |
---|
| 520 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi) |
---|
| 521 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) |
---|
| 522 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea) |
---|
| 523 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd) |
---|
| 524 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zsigS,zsig) |
---|
| 525 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam) |
---|
| 526 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe) |
---|
| 527 | |
---|
| 528 | endif ! of if (choix_1.eq.0) |
---|
| 529 | |
---|
| 530 | |
---|
| 531 | c======================================================================= |
---|
| 532 | c Lecture des fichiers (start ou start_archive) |
---|
| 533 | c======================================================================= |
---|
| 534 | |
---|
| 535 | if (choix_1.eq.0) then |
---|
| 536 | |
---|
| 537 | write(*,*) 'Reading file START_ARCHIVE' |
---|
[227] | 538 | CALL lect_start_archive(ngridmx,llm, |
---|
| 539 | & date,tsurf,tsoil,emis,q2, |
---|
[222] | 540 | & t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf, |
---|
| 541 | & surfith,nid, |
---|
| 542 | & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) |
---|
| 543 | write(*,*) "OK, read start_archive file" |
---|
| 544 | ! copy soil thermal inertia |
---|
| 545 | ithfi(:,:)=inertiedat(:,:) |
---|
| 546 | |
---|
| 547 | ierr= NF_CLOSE(nid) |
---|
| 548 | |
---|
| 549 | else if (choix_1.eq.1) then |
---|
| 550 | !do nothing, start and startfi have already been read |
---|
| 551 | else |
---|
| 552 | CALL exit(1) |
---|
| 553 | endif |
---|
| 554 | |
---|
| 555 | dtvr = daysec/FLOAT(day_step) |
---|
| 556 | dtphys = dtvr * FLOAT(iphysiq) |
---|
| 557 | |
---|
| 558 | c======================================================================= |
---|
| 559 | c |
---|
| 560 | c======================================================================= |
---|
| 561 | |
---|
| 562 | do ! infinite loop on list of changes |
---|
| 563 | |
---|
| 564 | write(*,*) |
---|
| 565 | write(*,*) |
---|
| 566 | write(*,*) 'List of possible changes :' |
---|
| 567 | write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
| 568 | write(*,*) |
---|
| 569 | write(*,*) 'flat : no topography ("aquaplanet")' |
---|
| 570 | write(*,*) 'nuketharsis : no Tharsis bulge' |
---|
| 571 | write(*,*) 'bilball : uniform albedo and thermal inertia' |
---|
| 572 | write(*,*) 'coldspole : cold subsurface and high albedo at S.pole' |
---|
| 573 | write(*,*) 'qname : change tracer name' |
---|
| 574 | write(*,*) 't=profile : read temperature profile in profile.in' |
---|
| 575 | write(*,*) 'q=0 : ALL tracer =zero' |
---|
| 576 | write(*,*) 'q=x : give a specific uniform value to one tracer' |
---|
| 577 | write(*,*) 'q=profile : specify a profile for a tracer' |
---|
| 578 | ! write(*,*) 'ini_q : tracers initialisation for chemistry, water an |
---|
| 579 | ! $d ice ' |
---|
| 580 | ! write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and |
---|
| 581 | ! $ice ' |
---|
| 582 | ! write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on |
---|
| 583 | ! $ly ' |
---|
| 584 | write(*,*) 'noglacier : Remove tropical H2O ice if |lat|<45' |
---|
| 585 | write(*,*) 'watercapn : H20 ice on permanent N polar cap ' |
---|
| 586 | write(*,*) 'watercaps : H20 ice on permanent S polar cap ' |
---|
| 587 | write(*,*) 'noacglac : H2O ice across Noachis Terra' |
---|
| 588 | write(*,*) 'oborealis : H2O ice across Vastitas Borealis' |
---|
| 589 | write(*,*) 'iceball : Thick ice layer all over surface' |
---|
| 590 | write(*,*) 'wetstart : start with a wet atmosphere' |
---|
| 591 | write(*,*) 'isotherm : Isothermal Temperatures, wind set to zero' |
---|
| 592 | write(*,*) 'radequi : Earth-like radiative equilibrium temperature |
---|
| 593 | $ profile (lat-alt) and winds set to zero' |
---|
| 594 | write(*,*) 'coldstart : Start X K above the CO2 frost point and |
---|
| 595 | $set wind to zero (assumes 100% CO2)' |
---|
| 596 | write(*,*) 'co2ice=0 : remove CO2 polar cap' |
---|
| 597 | write(*,*) 'ptot : change total pressure' |
---|
| 598 | write(*,*) 'emis : change surface emissivity' |
---|
| 599 | write(*,*) 'therm_ini_s : Set soil thermal inertia to reference sur |
---|
| 600 | &face values' |
---|
| 601 | write(*,*) 'slab_ocean_0 : initialisation of slab |
---|
| 602 | $ocean variables' |
---|
| 603 | |
---|
| 604 | write(*,*) |
---|
| 605 | write(*,*) 'Change to perform ?' |
---|
| 606 | write(*,*) ' (enter keyword or return to end)' |
---|
| 607 | write(*,*) |
---|
| 608 | |
---|
| 609 | read(*,fmt='(a20)') modif |
---|
| 610 | if (modif(1:1) .eq. ' ')then |
---|
| 611 | exit ! exit loop on changes |
---|
| 612 | endif |
---|
| 613 | |
---|
| 614 | write(*,*) |
---|
| 615 | write(*,*) trim(modif) , ' : ' |
---|
| 616 | |
---|
| 617 | c 'flat : no topography ("aquaplanet")' |
---|
| 618 | c ------------------------------------- |
---|
| 619 | if (trim(modif) .eq. 'flat') then |
---|
| 620 | c set topo to zero |
---|
| 621 | z_reel(1:iip1,1:jjp1)=0 |
---|
| 622 | phis(1:iip1,1:jjp1)=g*z_reel(1:iip1,1:jjp1) |
---|
| 623 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) |
---|
| 624 | write(*,*) 'topography set to zero.' |
---|
| 625 | write(*,*) 'WARNING : the subgrid topography parameters', |
---|
| 626 | & ' were not set to zero ! => set calllott to F' |
---|
| 627 | |
---|
| 628 | c Choice of surface pressure |
---|
| 629 | yes=' ' |
---|
| 630 | do while ((yes.ne.'y').and.(yes.ne.'n')) |
---|
| 631 | write(*,*) 'Do you wish to choose homogeneous surface', |
---|
| 632 | & 'pressure (y) or let newstart interpolate ', |
---|
| 633 | & ' the previous field (n)?' |
---|
| 634 | read(*,fmt='(a)') yes |
---|
| 635 | end do |
---|
| 636 | if (yes.eq.'y') then |
---|
| 637 | flagps0=.true. |
---|
| 638 | write(*,*) 'New value for ps (Pa) ?' |
---|
| 639 | 201 read(*,*,iostat=ierr) patm |
---|
| 640 | if(ierr.ne.0) goto 201 |
---|
| 641 | write(*,*) |
---|
| 642 | write(*,*) ' new ps everywhere (Pa) = ', patm |
---|
| 643 | write(*,*) |
---|
| 644 | do j=1,jjp1 |
---|
| 645 | do i=1,iip1 |
---|
| 646 | ps(i,j)=patm |
---|
| 647 | enddo |
---|
| 648 | enddo |
---|
| 649 | end if |
---|
| 650 | |
---|
| 651 | c 'nuketharsis : no tharsis bulge for Early Mars' |
---|
| 652 | c --------------------------------------------- |
---|
| 653 | else if (trim(modif) .eq. 'nuketharsis') then |
---|
| 654 | |
---|
| 655 | DO j=1,jjp1 |
---|
| 656 | DO i=1,iim |
---|
| 657 | ig=1+(j-2)*iim +i |
---|
| 658 | if(j.eq.1) ig=1 |
---|
| 659 | if(j.eq.jjp1) ig=ngridmx |
---|
| 660 | |
---|
| 661 | fact1=(((rlonv(i)*180./pi)+100)**2 + |
---|
| 662 | & (rlatu(j)*180./pi)**2)/65**2 |
---|
| 663 | fact2=exp( -fact1**2.5 ) |
---|
| 664 | |
---|
| 665 | phis(i,j) = phis(i,j) - (phis(i,j)+4000.*g)*fact2 |
---|
| 666 | |
---|
| 667 | ! if(phis(i,j).gt.2500.*g)then |
---|
| 668 | ! if(rlatu(j)*180./pi.gt.-80.)then ! avoid chopping south polar cap |
---|
| 669 | ! phis(i,j)=2500.*g |
---|
| 670 | ! endif |
---|
| 671 | ! endif |
---|
| 672 | |
---|
| 673 | ENDDO |
---|
| 674 | ENDDO |
---|
| 675 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) |
---|
| 676 | |
---|
| 677 | |
---|
| 678 | c bilball : uniform albedo, thermal inertia |
---|
| 679 | c ----------------------------------------- |
---|
| 680 | else if (trim(modif) .eq. 'bilball') then |
---|
| 681 | write(*,*) 'constante albedo and iner.therm:' |
---|
| 682 | write(*,*) 'New value for albedo (ex: 0.25) ?' |
---|
| 683 | 101 read(*,*,iostat=ierr) alb_bb |
---|
| 684 | if(ierr.ne.0) goto 101 |
---|
| 685 | write(*,*) |
---|
| 686 | write(*,*) ' uniform albedo (new value):',alb_bb |
---|
| 687 | write(*,*) |
---|
| 688 | |
---|
| 689 | write(*,*) 'New value for thermal inertia (eg: 247) ?' |
---|
| 690 | 102 read(*,*,iostat=ierr) ith_bb |
---|
| 691 | if(ierr.ne.0) goto 102 |
---|
| 692 | write(*,*) 'uniform thermal inertia (new value):',ith_bb |
---|
| 693 | DO j=1,jjp1 |
---|
| 694 | DO i=1,iip1 |
---|
| 695 | alb(i,j) = alb_bb ! albedo |
---|
| 696 | do isoil=1,nsoilmx |
---|
| 697 | ith(i,j,isoil) = ith_bb ! thermal inertia |
---|
| 698 | enddo |
---|
| 699 | END DO |
---|
| 700 | END DO |
---|
| 701 | ! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 702 | CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 703 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) |
---|
| 704 | |
---|
| 705 | c coldspole : sous-sol de la calotte sud toujours froid |
---|
| 706 | c ----------------------------------------------------- |
---|
| 707 | else if (trim(modif) .eq. 'coldspole') then |
---|
| 708 | write(*,*)'new value for the subsurface temperature', |
---|
| 709 | & ' beneath the permanent southern polar cap ? (eg: 141 K)' |
---|
| 710 | 103 read(*,*,iostat=ierr) tsud |
---|
| 711 | if(ierr.ne.0) goto 103 |
---|
| 712 | write(*,*) |
---|
| 713 | write(*,*) ' new value of the subsurface temperature:',tsud |
---|
| 714 | c nouvelle temperature sous la calotte permanente |
---|
| 715 | do l=2,nsoilmx |
---|
| 716 | tsoil(ngridmx,l) = tsud |
---|
| 717 | end do |
---|
| 718 | |
---|
| 719 | |
---|
| 720 | write(*,*)'new value for the albedo', |
---|
| 721 | & 'of the permanent southern polar cap ? (eg: 0.75)' |
---|
| 722 | 104 read(*,*,iostat=ierr) albsud |
---|
| 723 | if(ierr.ne.0) goto 104 |
---|
| 724 | write(*,*) |
---|
| 725 | |
---|
| 726 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 727 | c Option 1: only the albedo of the pole is modified : |
---|
| 728 | albfi(ngridmx)=albsud |
---|
| 729 | write(*,*) 'ig=',ngridmx,' albedo perennial cap ', |
---|
| 730 | & albfi(ngridmx) |
---|
| 731 | |
---|
| 732 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 733 | c Option 2 A haute resolution : coordonnee de la vrai calotte ~ |
---|
| 734 | c DO j=1,jjp1 |
---|
| 735 | c DO i=1,iip1 |
---|
| 736 | c ig=1+(j-2)*iim +i |
---|
| 737 | c if(j.eq.1) ig=1 |
---|
| 738 | c if(j.eq.jjp1) ig=ngridmx |
---|
| 739 | c if ((rlatu(j)*180./pi.lt.-84.).and. |
---|
| 740 | c & (rlatu(j)*180./pi.gt.-91.).and. |
---|
| 741 | c & (rlonv(i)*180./pi.gt.-91.).and. |
---|
| 742 | c & (rlonv(i)*180./pi.lt.0.)) then |
---|
| 743 | cc albedo de la calotte permanente fixe a albsud |
---|
| 744 | c alb(i,j)=albsud |
---|
| 745 | c write(*,*) 'lat=',rlatu(j)*180./pi, |
---|
| 746 | c & ' lon=',rlonv(i)*180./pi |
---|
| 747 | cc fin de la condition sur les limites de la calotte permanente |
---|
| 748 | c end if |
---|
| 749 | c ENDDO |
---|
| 750 | c ENDDO |
---|
| 751 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 752 | |
---|
| 753 | c CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) |
---|
| 754 | |
---|
| 755 | |
---|
| 756 | c ptot : Modification of the total pressure: ice + current atmosphere |
---|
| 757 | c ------------------------------------------------------------------- |
---|
| 758 | else if (trim(modif).eq.'ptot') then |
---|
| 759 | |
---|
| 760 | ! check if we have a co2_ice surface tracer: |
---|
| 761 | if (igcm_co2_ice.eq.0) then |
---|
| 762 | write(*,*) " No surface CO2 ice !" |
---|
| 763 | write(*,*) " only atmospheric pressure will be considered!" |
---|
| 764 | endif |
---|
| 765 | c calcul de la pression totale glace + atm actuelle |
---|
| 766 | patm=0. |
---|
| 767 | airetot=0. |
---|
| 768 | pcap=0. |
---|
| 769 | DO j=1,jjp1 |
---|
| 770 | DO i=1,iim |
---|
| 771 | ig=1+(j-2)*iim +i |
---|
| 772 | if(j.eq.1) ig=1 |
---|
| 773 | if(j.eq.jjp1) ig=ngridmx |
---|
| 774 | patm = patm + ps(i,j)*aire(i,j) |
---|
| 775 | airetot= airetot + aire(i,j) |
---|
| 776 | if (igcm_co2_ice.ne.0) then |
---|
| 777 | !pcap = pcap + aire(i,j)*co2ice(ig)*g |
---|
| 778 | pcap = pcap + aire(i,j)*qsurf(ig,igcm_co2_ice)*g |
---|
| 779 | endif |
---|
| 780 | ENDDO |
---|
| 781 | ENDDO |
---|
| 782 | ptoto = pcap + patm |
---|
| 783 | |
---|
| 784 | print*,'Current total pressure at surface (co2 ice + atm) ', |
---|
| 785 | & ptoto/airetot |
---|
| 786 | |
---|
| 787 | print*,'new value?' |
---|
| 788 | read(*,*) ptotn |
---|
| 789 | ptotn=ptotn*airetot |
---|
| 790 | patmn=ptotn-pcap |
---|
| 791 | print*,'ptoto,patm,ptotn,patmn' |
---|
| 792 | print*,ptoto,patm,ptotn,patmn |
---|
| 793 | print*,'Mult. factor for pressure (atm only)', patmn/patm |
---|
| 794 | do j=1,jjp1 |
---|
| 795 | do i=1,iip1 |
---|
| 796 | ps(i,j)=ps(i,j)*patmn/patm |
---|
| 797 | enddo |
---|
| 798 | enddo |
---|
| 799 | |
---|
| 800 | |
---|
| 801 | |
---|
| 802 | c Correction pour la conservation des traceurs |
---|
| 803 | yes=' ' |
---|
| 804 | do while ((yes.ne.'y').and.(yes.ne.'n')) |
---|
| 805 | write(*,*) 'Do you wish to conserve tracer total mass (y)', |
---|
| 806 | & ' or tracer mixing ratio (n) ?' |
---|
| 807 | read(*,fmt='(a)') yes |
---|
| 808 | end do |
---|
| 809 | |
---|
| 810 | if (yes.eq.'y') then |
---|
| 811 | write(*,*) 'OK : conservation of tracer total mass' |
---|
| 812 | DO iq =1, nqtot |
---|
| 813 | DO l=1,llm |
---|
| 814 | DO j=1,jjp1 |
---|
| 815 | DO i=1,iip1 |
---|
| 816 | q(i,j,l,iq)=q(i,j,l,iq)*patm/patmn |
---|
| 817 | ENDDO |
---|
| 818 | ENDDO |
---|
| 819 | ENDDO |
---|
| 820 | ENDDO |
---|
| 821 | else |
---|
| 822 | write(*,*) 'OK : conservation of tracer mixing ratio' |
---|
| 823 | end if |
---|
| 824 | |
---|
| 825 | c Correction pour la pression au niveau de la mer |
---|
| 826 | yes=' ' |
---|
| 827 | do while ((yes.ne.'y').and.(yes.ne.'n')) |
---|
| 828 | write(*,*) 'Do you wish fix pressure at sea level (y)', |
---|
| 829 | & ' or not (n) ?' |
---|
| 830 | read(*,fmt='(a)') yes |
---|
| 831 | end do |
---|
| 832 | |
---|
| 833 | if (yes.eq.'y') then |
---|
| 834 | write(*,*) 'Value?' |
---|
| 835 | read(*,*,iostat=ierr) psea |
---|
| 836 | DO i=1,iip1 |
---|
| 837 | DO j=1,jjp1 |
---|
| 838 | ps(i,j)=psea |
---|
| 839 | |
---|
| 840 | ENDDO |
---|
| 841 | ENDDO |
---|
| 842 | write(*,*) 'psea=',psea |
---|
| 843 | else |
---|
| 844 | write(*,*) 'no' |
---|
| 845 | end if |
---|
| 846 | |
---|
| 847 | |
---|
| 848 | c emis : change surface emissivity (added by RW) |
---|
| 849 | c ---------------------------------------------- |
---|
| 850 | else if (trim(modif).eq.'emis') then |
---|
| 851 | |
---|
| 852 | print*,'new value?' |
---|
| 853 | read(*,*) emisread |
---|
| 854 | |
---|
| 855 | do i=1,ngridmx |
---|
| 856 | emis(i)=emisread |
---|
| 857 | enddo |
---|
| 858 | |
---|
| 859 | c qname : change tracer name |
---|
| 860 | c -------------------------- |
---|
| 861 | else if (trim(modif).eq.'qname') then |
---|
| 862 | yes='y' |
---|
| 863 | do while (yes.eq.'y') |
---|
| 864 | write(*,*) 'Which tracer name do you want to change ?' |
---|
| 865 | do iq=1,nqtot |
---|
| 866 | write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq)) |
---|
| 867 | enddo |
---|
| 868 | write(*,'(a35,i3)') |
---|
| 869 | & '(enter tracer number; between 1 and ',nqtot |
---|
| 870 | write(*,*)' or any other value to quit this option)' |
---|
| 871 | read(*,*) iq |
---|
| 872 | if ((iq.ge.1).and.(iq.le.nqtot)) then |
---|
| 873 | write(*,*)'Change tracer name ',trim(tname(iq)),' to ?' |
---|
| 874 | read(*,*) txt |
---|
| 875 | tname(iq)=txt |
---|
| 876 | write(*,*)'Do you want to change another tracer name (y/n)?' |
---|
| 877 | read(*,'(a)') yes |
---|
| 878 | else |
---|
| 879 | ! inapropiate value of iq; quit this option |
---|
| 880 | yes='n' |
---|
| 881 | endif ! of if ((iq.ge.1).and.(iq.le.nqtot)) |
---|
| 882 | enddo ! of do while (yes.ne.'y') |
---|
| 883 | |
---|
| 884 | c q=0 : set tracers to zero |
---|
| 885 | c ------------------------- |
---|
| 886 | else if (trim(modif).eq.'q=0') then |
---|
| 887 | c mise a 0 des q (traceurs) |
---|
| 888 | write(*,*) 'Tracers set to 0 (1.E-30 in fact)' |
---|
| 889 | DO iq =1, nqtot |
---|
| 890 | DO l=1,llm |
---|
| 891 | DO j=1,jjp1 |
---|
| 892 | DO i=1,iip1 |
---|
| 893 | q(i,j,l,iq)=1.e-30 |
---|
| 894 | ENDDO |
---|
| 895 | ENDDO |
---|
| 896 | ENDDO |
---|
| 897 | ENDDO |
---|
| 898 | |
---|
| 899 | c set surface tracers to zero |
---|
| 900 | DO iq =1, nqtot |
---|
| 901 | DO ig=1,ngridmx |
---|
| 902 | qsurf(ig,iq)=0. |
---|
| 903 | ENDDO |
---|
| 904 | ENDDO |
---|
| 905 | |
---|
| 906 | c q=x : initialise tracer manually |
---|
| 907 | c -------------------------------- |
---|
| 908 | else if (trim(modif).eq.'q=x') then |
---|
| 909 | write(*,*) 'Which tracer do you want to modify ?' |
---|
| 910 | do iq=1,nqtot |
---|
| 911 | write(*,*)iq,' : ',trim(tname(iq)) |
---|
| 912 | enddo |
---|
| 913 | write(*,*) '(choose between 1 and ',nqtot,')' |
---|
| 914 | read(*,*) iq |
---|
| 915 | write(*,*)'mixing ratio of tracer ',trim(tname(iq)), |
---|
| 916 | & ' ? (kg/kg)' |
---|
| 917 | read(*,*) val |
---|
| 918 | DO l=1,llm |
---|
| 919 | DO j=1,jjp1 |
---|
| 920 | DO i=1,iip1 |
---|
| 921 | q(i,j,l,iq)=val |
---|
| 922 | ENDDO |
---|
| 923 | ENDDO |
---|
| 924 | ENDDO |
---|
| 925 | write(*,*) 'SURFACE value of tracer ',trim(tname(iq)), |
---|
| 926 | & ' ? (kg/m2)' |
---|
| 927 | read(*,*) val |
---|
| 928 | DO ig=1,ngridmx |
---|
| 929 | qsurf(ig,iq)=val |
---|
| 930 | ENDDO |
---|
| 931 | |
---|
| 932 | c t=profile : initialize temperature with a given profile |
---|
| 933 | c ------------------------------------------------------- |
---|
| 934 | else if (trim(modif) .eq. 't=profile') then |
---|
| 935 | write(*,*) 'Temperature profile from ASCII file' |
---|
| 936 | write(*,*) "'profile.in' e.g. 1D output" |
---|
| 937 | write(*,*) "(one value per line in file; starting with" |
---|
| 938 | write(*,*) "surface value, the 1st atmospheric layer" |
---|
| 939 | write(*,*) "followed by 2nd, etc. up to top of atmosphere)" |
---|
| 940 | txt="profile.in" |
---|
| 941 | open(33,file=trim(txt),status='old',form='formatted', |
---|
| 942 | & iostat=ierr) |
---|
| 943 | if (ierr.eq.0) then |
---|
| 944 | ! OK, found file 'profile_...', load the profile |
---|
| 945 | do l=1,llm+1 |
---|
| 946 | read(33,*,iostat=ierr) profile(l) |
---|
| 947 | write(*,*) profile(l) |
---|
| 948 | if (ierr.ne.0) then ! something went wrong |
---|
| 949 | exit ! quit loop |
---|
| 950 | endif |
---|
| 951 | enddo |
---|
| 952 | if (ierr.eq.0) then |
---|
| 953 | tsurf(1:ngridmx)=profile(1) |
---|
| 954 | tsoil(1:ngridmx,1:nsoilmx)=profile(1) |
---|
| 955 | do l=1,llm |
---|
| 956 | Tset(1:iip1,1:jjp1,l)=profile(l+1) |
---|
| 957 | flagtset=.true. |
---|
| 958 | enddo |
---|
| 959 | ucov(1:iip1,1:jjp1,1:llm)=0. |
---|
| 960 | vcov(1:iip1,1:jjm,1:llm)=0. |
---|
[227] | 961 | q2(1:ngridmx,1:llm+1)=0. |
---|
[222] | 962 | else |
---|
| 963 | write(*,*)'problem reading file ',trim(txt),' !' |
---|
| 964 | write(*,*)'No modifications to temperature' |
---|
| 965 | endif |
---|
| 966 | else |
---|
| 967 | write(*,*)'Could not find file ',trim(txt),' !' |
---|
| 968 | write(*,*)'No modifications to temperature' |
---|
| 969 | endif |
---|
| 970 | |
---|
| 971 | c q=profile : initialize tracer with a given profile |
---|
| 972 | c -------------------------------------------------- |
---|
| 973 | else if (trim(modif) .eq. 'q=profile') then |
---|
| 974 | write(*,*) 'Tracer profile will be sought in ASCII file' |
---|
| 975 | write(*,*) "'profile_tracer' where 'tracer' is tracer name" |
---|
| 976 | write(*,*) "(one value per line in file; starting with" |
---|
| 977 | write(*,*) "surface value, the 1st atmospheric layer" |
---|
| 978 | write(*,*) "followed by 2nd, etc. up to top of atmosphere)" |
---|
| 979 | write(*,*) 'Which tracer do you want to set?' |
---|
| 980 | do iq=1,nqtot |
---|
| 981 | write(*,*)iq,' : ',trim(tname(iq)) |
---|
| 982 | enddo |
---|
| 983 | write(*,*) '(choose between 1 and ',nqtot,')' |
---|
| 984 | read(*,*) iq |
---|
| 985 | if ((iq.lt.1).or.(iq.gt.nqtot)) then |
---|
| 986 | ! wrong value for iq, go back to menu |
---|
| 987 | write(*,*) "wrong input value:",iq |
---|
| 988 | cycle |
---|
| 989 | endif |
---|
| 990 | ! look for input file 'profile_tracer' |
---|
| 991 | txt="profile_"//trim(tname(iq)) |
---|
| 992 | open(41,file=trim(txt),status='old',form='formatted', |
---|
| 993 | & iostat=ierr) |
---|
| 994 | if (ierr.eq.0) then |
---|
| 995 | ! OK, found file 'profile_...', load the profile |
---|
| 996 | do l=1,llm+1 |
---|
| 997 | read(41,*,iostat=ierr) profile(l) |
---|
| 998 | if (ierr.ne.0) then ! something went wrong |
---|
| 999 | exit ! quit loop |
---|
| 1000 | endif |
---|
| 1001 | enddo |
---|
| 1002 | if (ierr.eq.0) then |
---|
| 1003 | ! initialize tracer values |
---|
| 1004 | qsurf(:,iq)=profile(1) |
---|
| 1005 | do l=1,llm |
---|
| 1006 | q(:,:,l,iq)=profile(l+1) |
---|
| 1007 | enddo |
---|
| 1008 | write(*,*)'OK, tracer ',trim(tname(iq)),' initialized ' |
---|
| 1009 | & ,'using values from file ',trim(txt) |
---|
| 1010 | else |
---|
| 1011 | write(*,*)'problem reading file ',trim(txt),' !' |
---|
| 1012 | write(*,*)'No modifications to tracer ',trim(tname(iq)) |
---|
| 1013 | endif |
---|
| 1014 | else |
---|
| 1015 | write(*,*)'Could not find file ',trim(txt),' !' |
---|
| 1016 | write(*,*)'No modifications to tracer ',trim(tname(iq)) |
---|
| 1017 | endif |
---|
| 1018 | |
---|
| 1019 | |
---|
| 1020 | c wetstart : wet atmosphere with a north to south gradient |
---|
| 1021 | c -------------------------------------------------------- |
---|
| 1022 | else if (trim(modif) .eq. 'wetstart') then |
---|
| 1023 | ! check that there is indeed a water vapor tracer |
---|
| 1024 | if (igcm_h2o_vap.eq.0) then |
---|
| 1025 | write(*,*) "No water vapour tracer! Can't use this option" |
---|
| 1026 | stop |
---|
| 1027 | endif |
---|
| 1028 | DO l=1,llm |
---|
| 1029 | DO j=1,jjp1 |
---|
| 1030 | DO i=1,iip1 |
---|
| 1031 | q(i,j,l,igcm_h2o_vap)=150.e-6 * (rlatu(j)+pi/2.) / pi |
---|
| 1032 | ENDDO |
---|
| 1033 | ENDDO |
---|
| 1034 | ENDDO |
---|
| 1035 | |
---|
| 1036 | write(*,*) 'Water mass mixing ratio at north pole=' |
---|
| 1037 | * ,q(1,1,1,igcm_h2o_vap) |
---|
| 1038 | write(*,*) '---------------------------south pole=' |
---|
| 1039 | * ,q(1,jjp1,1,igcm_h2o_vap) |
---|
| 1040 | |
---|
| 1041 | c noglacier : remove tropical water ice (to initialize high res sim) |
---|
| 1042 | c -------------------------------------------------- |
---|
| 1043 | else if (trim(modif) .eq. 'noglacier') then |
---|
| 1044 | if (igcm_h2o_ice.eq.0) then |
---|
| 1045 | write(*,*) "No water ice tracer! Can't use this option" |
---|
| 1046 | stop |
---|
| 1047 | endif |
---|
| 1048 | do ig=1,ngridmx |
---|
| 1049 | j=(ig-2)/iim +2 |
---|
| 1050 | if(ig.eq.1) j=1 |
---|
| 1051 | write(*,*) 'OK: remove surface ice for |lat|<45' |
---|
| 1052 | if (abs(rlatu(j)*180./pi).lt.45.) then |
---|
| 1053 | qsurf(ig,igcm_h2o_ice)=0. |
---|
| 1054 | end if |
---|
| 1055 | end do |
---|
| 1056 | |
---|
| 1057 | |
---|
| 1058 | c watercapn : H20 ice on permanent northern cap |
---|
| 1059 | c -------------------------------------------------- |
---|
| 1060 | else if (trim(modif) .eq. 'watercapn') then |
---|
| 1061 | if (igcm_h2o_ice.eq.0) then |
---|
| 1062 | write(*,*) "No water ice tracer! Can't use this option" |
---|
| 1063 | stop |
---|
| 1064 | endif |
---|
| 1065 | |
---|
| 1066 | DO j=1,jjp1 |
---|
| 1067 | DO i=1,iim |
---|
| 1068 | ig=1+(j-2)*iim +i |
---|
| 1069 | if(j.eq.1) ig=1 |
---|
| 1070 | if(j.eq.jjp1) ig=ngridmx |
---|
| 1071 | |
---|
| 1072 | if (rlatu(j)*180./pi.gt.80.) then |
---|
| 1073 | qsurf(ig,igcm_h2o_ice)=3.4e3 |
---|
| 1074 | !do isoil=1,nsoilmx |
---|
| 1075 | ! ith(i,j,isoil) = 18000. ! thermal inertia |
---|
| 1076 | !enddo |
---|
| 1077 | write(*,*)' ==> Ice mesh North boundary (deg)= ', |
---|
| 1078 | & rlatv(j-1)*180./pi |
---|
| 1079 | end if |
---|
| 1080 | ENDDO |
---|
| 1081 | ENDDO |
---|
| 1082 | CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 1083 | |
---|
| 1084 | c$$$ do ig=1,ngridmx |
---|
| 1085 | c$$$ j=(ig-2)/iim +2 |
---|
| 1086 | c$$$ if(ig.eq.1) j=1 |
---|
| 1087 | c$$$ if (rlatu(j)*180./pi.gt.80.) then |
---|
| 1088 | c$$$ |
---|
| 1089 | c$$$ qsurf(ig,igcm_h2o_ice)=1.e5 |
---|
| 1090 | c$$$ qsurf(ig,igcm_h2o_vap)=0.0!1.e5 |
---|
| 1091 | c$$$ |
---|
| 1092 | c$$$ write(*,*) 'ig=',ig,' H2O ice mass (kg/m2)= ', |
---|
| 1093 | c$$$ & qsurf(ig,igcm_h2o_ice) |
---|
| 1094 | c$$$ |
---|
| 1095 | c$$$ write(*,*)' ==> Ice mesh South boundary (deg)= ', |
---|
| 1096 | c$$$ & rlatv(j)*180./pi |
---|
| 1097 | c$$$ end if |
---|
| 1098 | c$$$ enddo |
---|
| 1099 | |
---|
| 1100 | c watercaps : H20 ice on permanent southern cap |
---|
| 1101 | c ------------------------------------------------- |
---|
| 1102 | else if (trim(modif) .eq. 'watercaps') then |
---|
| 1103 | if (igcm_h2o_ice.eq.0) then |
---|
| 1104 | write(*,*) "No water ice tracer! Can't use this option" |
---|
| 1105 | stop |
---|
| 1106 | endif |
---|
| 1107 | |
---|
| 1108 | DO j=1,jjp1 |
---|
| 1109 | DO i=1,iim |
---|
| 1110 | ig=1+(j-2)*iim +i |
---|
| 1111 | if(j.eq.1) ig=1 |
---|
| 1112 | if(j.eq.jjp1) ig=ngridmx |
---|
| 1113 | |
---|
| 1114 | if (rlatu(j)*180./pi.lt.-80.) then |
---|
| 1115 | qsurf(ig,igcm_h2o_ice)=3.4e3 |
---|
| 1116 | !do isoil=1,nsoilmx |
---|
| 1117 | ! ith(i,j,isoil) = 18000. ! thermal inertia |
---|
| 1118 | !enddo |
---|
| 1119 | write(*,*)' ==> Ice mesh South boundary (deg)= ', |
---|
| 1120 | & rlatv(j-1)*180./pi |
---|
| 1121 | end if |
---|
| 1122 | ENDDO |
---|
| 1123 | ENDDO |
---|
| 1124 | CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 1125 | |
---|
| 1126 | c$$$ do ig=1,ngridmx |
---|
| 1127 | c$$$ j=(ig-2)/iim +2 |
---|
| 1128 | c$$$ if(ig.eq.1) j=1 |
---|
| 1129 | c$$$ if (rlatu(j)*180./pi.lt.-80.) then |
---|
| 1130 | c$$$ qsurf(ig,igcm_h2o_ice)=1.e5 |
---|
| 1131 | c$$$ qsurf(ig,igcm_h2o_vap)=0.0 !1.e5 |
---|
| 1132 | c$$$ |
---|
| 1133 | c$$$ write(*,*) 'ig=',ig,' H2O ice mass (kg/m2)= ', |
---|
| 1134 | c$$$ & qsurf(ig,igcm_h2o_ice) |
---|
| 1135 | c$$$ write(*,*)' ==> Ice mesh North boundary (deg)= ', |
---|
| 1136 | c$$$ & rlatv(j-1)*180./pi |
---|
| 1137 | c$$$ end if |
---|
| 1138 | c$$$ enddo |
---|
| 1139 | |
---|
| 1140 | |
---|
| 1141 | c noacglac : H2O ice across highest terrain |
---|
| 1142 | c -------------------------------------------- |
---|
| 1143 | else if (trim(modif) .eq. 'noacglac') then |
---|
| 1144 | if (igcm_h2o_ice.eq.0) then |
---|
| 1145 | write(*,*) "No water ice tracer! Can't use this option" |
---|
| 1146 | stop |
---|
| 1147 | endif |
---|
| 1148 | DO j=1,jjp1 |
---|
| 1149 | DO i=1,iim |
---|
| 1150 | ig=1+(j-2)*iim +i |
---|
| 1151 | if(j.eq.1) ig=1 |
---|
| 1152 | if(j.eq.jjp1) ig=ngridmx |
---|
| 1153 | |
---|
| 1154 | if(phis(i,j).gt.1000.*g)then |
---|
| 1155 | alb(i,j) = 0.5 ! snow value |
---|
| 1156 | do isoil=1,nsoilmx |
---|
| 1157 | ith(i,j,isoil) = 18000. ! thermal inertia |
---|
| 1158 | ! this leads to rnat set to 'ocean' in physiq.F90 |
---|
| 1159 | ! actually no, because it is soil not surface |
---|
| 1160 | enddo |
---|
| 1161 | endif |
---|
| 1162 | ENDDO |
---|
| 1163 | ENDDO |
---|
| 1164 | CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 1165 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) |
---|
| 1166 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) |
---|
| 1167 | |
---|
| 1168 | |
---|
| 1169 | |
---|
| 1170 | c oborealis : H2O oceans across Vastitas Borealis |
---|
| 1171 | c ----------------------------------------------- |
---|
| 1172 | else if (trim(modif) .eq. 'oborealis') then |
---|
| 1173 | if (igcm_h2o_ice.eq.0) then |
---|
| 1174 | write(*,*) "No water ice tracer! Can't use this option" |
---|
| 1175 | stop |
---|
| 1176 | endif |
---|
| 1177 | DO j=1,jjp1 |
---|
| 1178 | DO i=1,iim |
---|
| 1179 | ig=1+(j-2)*iim +i |
---|
| 1180 | if(j.eq.1) ig=1 |
---|
| 1181 | if(j.eq.jjp1) ig=ngridmx |
---|
| 1182 | |
---|
| 1183 | if(phis(i,j).lt.-4000.*g)then |
---|
| 1184 | ! if( (phis(i,j).lt.-4000.*g) |
---|
| 1185 | ! & .and. (rlatu(j)*180./pi.lt.0.) )then ! south hemisphere only |
---|
| 1186 | ! phis(i,j)=-8000.0*g ! proper ocean |
---|
| 1187 | phis(i,j)=-4000.0*g ! scanty ocean |
---|
| 1188 | |
---|
| 1189 | alb(i,j) = 0.07 ! oceanic value |
---|
| 1190 | do isoil=1,nsoilmx |
---|
| 1191 | ith(i,j,isoil) = 18000. ! thermal inertia |
---|
| 1192 | ! this leads to rnat set to 'ocean' in physiq.F90 |
---|
| 1193 | ! actually no, because it is soil not surface |
---|
| 1194 | enddo |
---|
| 1195 | endif |
---|
| 1196 | ENDDO |
---|
| 1197 | ENDDO |
---|
| 1198 | CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 1199 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) |
---|
| 1200 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) |
---|
| 1201 | |
---|
| 1202 | c iborealis : H2O ice in Northern plains |
---|
| 1203 | c -------------------------------------- |
---|
| 1204 | else if (trim(modif) .eq. 'iborealis') then |
---|
| 1205 | if (igcm_h2o_ice.eq.0) then |
---|
| 1206 | write(*,*) "No water ice tracer! Can't use this option" |
---|
| 1207 | stop |
---|
| 1208 | endif |
---|
| 1209 | DO j=1,jjp1 |
---|
| 1210 | DO i=1,iim |
---|
| 1211 | ig=1+(j-2)*iim +i |
---|
| 1212 | if(j.eq.1) ig=1 |
---|
| 1213 | if(j.eq.jjp1) ig=ngridmx |
---|
| 1214 | |
---|
| 1215 | if(phis(i,j).lt.-4000.*g)then |
---|
| 1216 | !qsurf(ig,igcm_h2o_ice)=1.e3 |
---|
| 1217 | qsurf(ig,igcm_h2o_ice)=241.4 ! to make total 33 kg m^-2 |
---|
| 1218 | endif |
---|
| 1219 | ENDDO |
---|
| 1220 | ENDDO |
---|
| 1221 | CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 1222 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) |
---|
| 1223 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) |
---|
| 1224 | |
---|
| 1225 | |
---|
| 1226 | c oceanball : H2O liquid everywhere |
---|
| 1227 | c ---------------------------- |
---|
| 1228 | else if (trim(modif) .eq. 'oceanball') then |
---|
| 1229 | if (igcm_h2o_ice.eq.0) then |
---|
| 1230 | write(*,*) "No water ice tracer! Can't use this option" |
---|
| 1231 | stop |
---|
| 1232 | endif |
---|
| 1233 | DO j=1,jjp1 |
---|
| 1234 | DO i=1,iim |
---|
| 1235 | ig=1+(j-2)*iim +i |
---|
| 1236 | if(j.eq.1) ig=1 |
---|
| 1237 | if(j.eq.jjp1) ig=ngridmx |
---|
| 1238 | |
---|
| 1239 | qsurf(ig,igcm_h2o_vap)=0.0 ! for ocean, this is infinite source |
---|
| 1240 | qsurf(ig,igcm_h2o_ice)=0.0 |
---|
| 1241 | alb(i,j) = 0.07 ! ocean value |
---|
| 1242 | |
---|
| 1243 | do isoil=1,nsoilmx |
---|
| 1244 | ith(i,j,isoil) = 18000. ! thermal inertia |
---|
| 1245 | !ith(i,j,isoil) = 50. ! extremely low for test |
---|
| 1246 | ! this leads to rnat set to 'ocean' in physiq.F90 |
---|
| 1247 | enddo |
---|
| 1248 | |
---|
| 1249 | ENDDO |
---|
| 1250 | ENDDO |
---|
| 1251 | CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 1252 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) |
---|
| 1253 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi) |
---|
| 1254 | |
---|
| 1255 | c iceball : H2O ice everywhere |
---|
| 1256 | c ---------------------------- |
---|
| 1257 | else if (trim(modif) .eq. 'iceball') then |
---|
| 1258 | if (igcm_h2o_ice.eq.0) then |
---|
| 1259 | write(*,*) "No water ice tracer! Can't use this option" |
---|
| 1260 | stop |
---|
| 1261 | endif |
---|
| 1262 | DO j=1,jjp1 |
---|
| 1263 | DO i=1,iim |
---|
| 1264 | ig=1+(j-2)*iim +i |
---|
| 1265 | if(j.eq.1) ig=1 |
---|
| 1266 | if(j.eq.jjp1) ig=ngridmx |
---|
| 1267 | |
---|
| 1268 | qsurf(ig,igcm_h2o_vap)=-50. ! for ocean, this is infinite source |
---|
| 1269 | qsurf(ig,igcm_h2o_ice)=50. ! == to 0.05 m of oceanic ice |
---|
| 1270 | alb(i,j) = 0.6 ! ice albedo value |
---|
| 1271 | |
---|
| 1272 | do isoil=1,nsoilmx |
---|
| 1273 | !ith(i,j,isoil) = 18000. ! thermal inertia |
---|
| 1274 | ! this leads to rnat set to 'ocean' in physiq.F90 |
---|
| 1275 | enddo |
---|
| 1276 | |
---|
| 1277 | ENDDO |
---|
| 1278 | ENDDO |
---|
| 1279 | CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi) |
---|
| 1280 | CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi) |
---|
| 1281 | |
---|
| 1282 | c isotherm : Isothermal temperatures and no winds |
---|
| 1283 | c ----------------------------------------------- |
---|
| 1284 | else if (trim(modif) .eq. 'isotherm') then |
---|
| 1285 | |
---|
| 1286 | write(*,*)'Isothermal temperature of the atmosphere, |
---|
| 1287 | & surface and subsurface' |
---|
| 1288 | write(*,*) 'Value of this temperature ? :' |
---|
| 1289 | 203 read(*,*,iostat=ierr) Tiso |
---|
| 1290 | if(ierr.ne.0) goto 203 |
---|
| 1291 | |
---|
| 1292 | tsurf(1:ngridmx)=Tiso |
---|
| 1293 | |
---|
| 1294 | tsoil(1:ngridmx,1:nsoilmx)=Tiso |
---|
| 1295 | |
---|
| 1296 | Tset(1:iip1,1:jjp1,1:llm)=Tiso |
---|
| 1297 | flagtset=.true. |
---|
| 1298 | |
---|
| 1299 | ucov(1:iip1,1:jjp1,1:llm)=0 |
---|
| 1300 | vcov(1:iip1,1:jjm,1:llm)=0 |
---|
[227] | 1301 | q2(1:ngridmx,1:llm+1)=0 |
---|
[222] | 1302 | |
---|
| 1303 | c radequi : Radiative equilibrium profile of temperatures and no winds |
---|
| 1304 | c -------------------------------------------------------------------- |
---|
| 1305 | else if (trim(modif) .eq. 'radequi') then |
---|
| 1306 | |
---|
| 1307 | write(*,*)'radiative equilibrium temperature profile' |
---|
| 1308 | |
---|
| 1309 | do ig=1, ngridmx |
---|
| 1310 | teque= 335.0-60.0*sin(latfi(ig))*sin(latfi(ig))- |
---|
| 1311 | & 10.0*cos(latfi(ig))*cos(latfi(ig)) |
---|
| 1312 | tsurf(ig) = MAX(220.0,teque) |
---|
| 1313 | end do |
---|
| 1314 | do l=2,nsoilmx |
---|
| 1315 | do ig=1, ngridmx |
---|
| 1316 | tsoil(ig,l) = tsurf(ig) |
---|
| 1317 | end do |
---|
| 1318 | end do |
---|
| 1319 | DO j=1,jjp1 |
---|
| 1320 | DO i=1,iim |
---|
| 1321 | Do l=1,llm |
---|
| 1322 | teque=335.-60.0*sin(rlatu(j))*sin(rlatu(j))- |
---|
| 1323 | & 10.0*cos(rlatu(j))*cos(rlatu(j)) |
---|
| 1324 | Tset(i,j,l)=MAX(220.0,teque) |
---|
| 1325 | end do |
---|
| 1326 | end do |
---|
| 1327 | end do |
---|
| 1328 | flagtset=.true. |
---|
| 1329 | ucov(1:iip1,1:jjp1,1:llm)=0 |
---|
| 1330 | vcov(1:iip1,1:jjm,1:llm)=0 |
---|
[227] | 1331 | q2(1:ngridmx,1:llm+1)=0 |
---|
[222] | 1332 | |
---|
| 1333 | c coldstart : T set 1K above CO2 frost point and no winds |
---|
| 1334 | c ------------------------------------------------ |
---|
| 1335 | else if (trim(modif) .eq. 'coldstart') then |
---|
| 1336 | |
---|
| 1337 | write(*,*)'set temperature of the atmosphere,' |
---|
| 1338 | &,'surface and subsurface how many degrees above CO2 frost point?' |
---|
| 1339 | 204 read(*,*,iostat=ierr) Tabove |
---|
| 1340 | if(ierr.ne.0) goto 204 |
---|
| 1341 | |
---|
| 1342 | DO j=1,jjp1 |
---|
| 1343 | DO i=1,iim |
---|
| 1344 | ig=1+(j-2)*iim +i |
---|
| 1345 | if(j.eq.1) ig=1 |
---|
| 1346 | if(j.eq.jjp1) ig=ngridmx |
---|
| 1347 | tsurf(ig) = (-3167.8)/(log(.01*ps(i,j))-23.23)+Tabove |
---|
| 1348 | END DO |
---|
| 1349 | END DO |
---|
| 1350 | do l=1,nsoilmx |
---|
| 1351 | do ig=1, ngridmx |
---|
| 1352 | tsoil(ig,l) = tsurf(ig) |
---|
| 1353 | end do |
---|
| 1354 | end do |
---|
| 1355 | DO j=1,jjp1 |
---|
| 1356 | DO i=1,iim |
---|
| 1357 | Do l=1,llm |
---|
| 1358 | pp = aps(l) +bps(l)*ps(i,j) |
---|
| 1359 | Tset(i,j,l)=(-3167.8)/(log(.01*pp)-23.23)+Tabove |
---|
| 1360 | end do |
---|
| 1361 | end do |
---|
| 1362 | end do |
---|
| 1363 | |
---|
| 1364 | flagtset=.true. |
---|
| 1365 | ucov(1:iip1,1:jjp1,1:llm)=0 |
---|
| 1366 | vcov(1:iip1,1:jjm,1:llm)=0 |
---|
[227] | 1367 | q2(1:ngridmx,1:llm+1)=0 |
---|
[222] | 1368 | |
---|
| 1369 | |
---|
| 1370 | c co2ice=0 : remove CO2 polar ice caps' |
---|
| 1371 | c ------------------------------------------------ |
---|
| 1372 | else if (trim(modif) .eq. 'co2ice=0') then |
---|
| 1373 | ! check that there is indeed a co2_ice tracer ... |
---|
| 1374 | if (igcm_co2_ice.ne.0) then |
---|
| 1375 | do ig=1,ngridmx |
---|
| 1376 | !co2ice(ig)=0 |
---|
| 1377 | qsurf(ig,igcm_co2_ice)=0 |
---|
| 1378 | emis(ig)=emis(ngridmx/2) |
---|
| 1379 | end do |
---|
| 1380 | else |
---|
| 1381 | write(*,*) "Can't remove CO2 ice!! (no co2_ice tracer)" |
---|
| 1382 | endif |
---|
| 1383 | |
---|
| 1384 | ! therm_ini_s: (re)-set soil thermal inertia to reference surface values |
---|
| 1385 | ! ---------------------------------------------------------------------- |
---|
| 1386 | |
---|
| 1387 | else if (trim(modif) .eq. 'therm_ini_s') then |
---|
| 1388 | ! write(*,*)"surfithfi(1):",surfithfi(1) |
---|
| 1389 | do isoil=1,nsoilmx |
---|
| 1390 | inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx) |
---|
| 1391 | enddo |
---|
| 1392 | write(*,*)'OK: Soil thermal inertia has been reset to referenc |
---|
| 1393 | &e surface values' |
---|
| 1394 | ! write(*,*)"inertiedat(1,1):",inertiedat(1,1) |
---|
| 1395 | ithfi(:,:)=inertiedat(:,:) |
---|
| 1396 | ! recast ithfi() onto ith() |
---|
| 1397 | call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith) |
---|
| 1398 | ! Check: |
---|
| 1399 | ! do i=1,iip1 |
---|
| 1400 | ! do j=1,jjp1 |
---|
| 1401 | ! do isoil=1,nsoilmx |
---|
| 1402 | ! write(77,*) i,j,isoil," ",ith(i,j,isoil) |
---|
| 1403 | ! enddo |
---|
| 1404 | ! enddo |
---|
| 1405 | ! enddo |
---|
| 1406 | |
---|
| 1407 | |
---|
| 1408 | |
---|
| 1409 | c slab_ocean_initialisation |
---|
| 1410 | c ------------------------------------------------ |
---|
| 1411 | else if (trim(modif) .eq. 'slab_ocean_0') then |
---|
| 1412 | write(*,*)'OK: initialisation of slab ocean' |
---|
| 1413 | |
---|
| 1414 | DO ig=1, ngridmx |
---|
| 1415 | rnat(ig)=1. |
---|
| 1416 | tslab(ig,1)=0. |
---|
| 1417 | tslab(ig,2)=0. |
---|
| 1418 | tsea_ice(ig)=0. |
---|
| 1419 | sea_ice(ig)=0. |
---|
| 1420 | pctsrf_sic(ig)=0. |
---|
| 1421 | |
---|
| 1422 | if(ithfi(ig,1).GT.10000.)then |
---|
| 1423 | rnat(ig)=0. |
---|
| 1424 | phisfi(ig)=0. |
---|
| 1425 | tsea_ice(ig)=273.15-1.8 |
---|
| 1426 | tslab(ig,1)=tsurf(ig) |
---|
| 1427 | tslab(ig,2)=tsurf(ig)!*2./3.+(273.15-1.8)/3. |
---|
| 1428 | endif |
---|
| 1429 | |
---|
| 1430 | ENDDO |
---|
| 1431 | CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,phisfi,phis) |
---|
| 1432 | |
---|
| 1433 | |
---|
| 1434 | |
---|
| 1435 | else |
---|
| 1436 | write(*,*) ' Unknown (misspelled?) option!!!' |
---|
| 1437 | end if ! of if (trim(modif) .eq. '...') elseif ... |
---|
| 1438 | |
---|
| 1439 | |
---|
| 1440 | |
---|
| 1441 | enddo ! of do ! infinite loop on liste of changes |
---|
| 1442 | |
---|
| 1443 | 999 continue |
---|
| 1444 | |
---|
| 1445 | |
---|
| 1446 | c======================================================================= |
---|
| 1447 | c Initialisation for cloud fraction and oceanic ice (added by BC 2010) |
---|
| 1448 | c======================================================================= |
---|
| 1449 | DO ig=1, ngridmx |
---|
| 1450 | totalfrac(ig)=0.5 |
---|
| 1451 | DO l=1,llm |
---|
| 1452 | cloudfrac(ig,l)=0.5 |
---|
| 1453 | ENDDO |
---|
| 1454 | ! Ehouarn, march 2012: also add some initialisation for hice |
---|
| 1455 | hice(ig)=0.0 |
---|
| 1456 | ENDDO |
---|
| 1457 | |
---|
| 1458 | c======================================================== |
---|
| 1459 | |
---|
| 1460 | ! DO ig=1,ngridmx |
---|
| 1461 | ! IF(tsurf(ig) .LT. 273.16-1.8) THEN |
---|
| 1462 | ! hice(ig)=(273.16-1.8-tsurf(ig))/(273.16-1.8-240)*1. |
---|
| 1463 | ! hice(ig)=min(hice(ig),1.0) |
---|
| 1464 | ! ENDIF |
---|
| 1465 | ! ENDDO |
---|
| 1466 | |
---|
| 1467 | |
---|
| 1468 | |
---|
| 1469 | |
---|
| 1470 | c======================================================================= |
---|
| 1471 | c Correct pressure on the new grid (menu 0) |
---|
| 1472 | c======================================================================= |
---|
| 1473 | |
---|
| 1474 | |
---|
| 1475 | if ((choix_1.eq.0).and.(.not.flagps0)) then |
---|
| 1476 | r = 1000.*8.31/mugaz |
---|
| 1477 | |
---|
| 1478 | do j=1,jjp1 |
---|
| 1479 | do i=1,iip1 |
---|
| 1480 | ps(i,j) = ps(i,j) * |
---|
| 1481 | . exp((phisold_newgrid(i,j)-phis(i,j)) / |
---|
| 1482 | . (t(i,j,1) * r)) |
---|
| 1483 | end do |
---|
| 1484 | end do |
---|
| 1485 | |
---|
| 1486 | c periodicite de ps en longitude |
---|
| 1487 | do j=1,jjp1 |
---|
| 1488 | ps(1,j) = ps(iip1,j) |
---|
| 1489 | end do |
---|
| 1490 | end if |
---|
| 1491 | |
---|
| 1492 | |
---|
| 1493 | c======================================================================= |
---|
| 1494 | c======================================================================= |
---|
| 1495 | |
---|
| 1496 | c======================================================================= |
---|
| 1497 | c Initialisation de la physique / ecriture de newstartfi : |
---|
| 1498 | c======================================================================= |
---|
| 1499 | |
---|
| 1500 | |
---|
| 1501 | CALL inifilr |
---|
| 1502 | CALL pression(ip1jmp1, ap, bp, ps, p3d) |
---|
| 1503 | |
---|
| 1504 | c----------------------------------------------------------------------- |
---|
| 1505 | c Initialisation pks: |
---|
| 1506 | c----------------------------------------------------------------------- |
---|
| 1507 | |
---|
| 1508 | CALL exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf) |
---|
| 1509 | ! Calcul de la temperature potentielle teta |
---|
| 1510 | |
---|
| 1511 | if (flagtset) then |
---|
| 1512 | DO l=1,llm |
---|
| 1513 | DO j=1,jjp1 |
---|
| 1514 | DO i=1,iim |
---|
| 1515 | teta(i,j,l) = Tset(i,j,l) * cpp/pk(i,j,l) |
---|
| 1516 | ENDDO |
---|
| 1517 | teta (iip1,j,l)= teta (1,j,l) |
---|
| 1518 | ENDDO |
---|
| 1519 | ENDDO |
---|
| 1520 | else if (choix_1.eq.0) then |
---|
| 1521 | DO l=1,llm |
---|
| 1522 | DO j=1,jjp1 |
---|
| 1523 | DO i=1,iim |
---|
| 1524 | teta(i,j,l) = t(i,j,l) * cpp/pk(i,j,l) |
---|
| 1525 | ENDDO |
---|
| 1526 | teta (iip1,j,l)= teta (1,j,l) |
---|
| 1527 | ENDDO |
---|
| 1528 | ENDDO |
---|
| 1529 | endif |
---|
| 1530 | |
---|
| 1531 | C Calcul intermediaire |
---|
| 1532 | c |
---|
| 1533 | if (choix_1.eq.0) then |
---|
| 1534 | CALL massdair( p3d, masse ) |
---|
| 1535 | c |
---|
| 1536 | print *,' ALPHAX ',alphax |
---|
| 1537 | |
---|
| 1538 | DO l = 1, llm |
---|
| 1539 | DO i = 1, iim |
---|
| 1540 | xppn(i) = aire( i, 1 ) * masse( i , 1 , l ) |
---|
| 1541 | xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l ) |
---|
| 1542 | ENDDO |
---|
| 1543 | xpn = SUM(xppn)/apoln |
---|
| 1544 | xps = SUM(xpps)/apols |
---|
| 1545 | DO i = 1, iip1 |
---|
| 1546 | masse( i , 1 , l ) = xpn |
---|
| 1547 | masse( i , jjp1 , l ) = xps |
---|
| 1548 | ENDDO |
---|
| 1549 | ENDDO |
---|
| 1550 | endif |
---|
| 1551 | phis(iip1,:) = phis(1,:) |
---|
| 1552 | |
---|
| 1553 | itau=0 |
---|
| 1554 | if (choix_1.eq.0) then |
---|
| 1555 | day_ini=int(date) |
---|
| 1556 | endif |
---|
| 1557 | c |
---|
| 1558 | CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) |
---|
| 1559 | |
---|
| 1560 | CALL caldyn0( itau,ucov,vcov,teta,ps,masse,pk,phis , |
---|
| 1561 | * phi,w, pbaru,pbarv,day_ini+time ) |
---|
| 1562 | |
---|
| 1563 | |
---|
| 1564 | CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqtot) |
---|
| 1565 | CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqtot,masse,ps) |
---|
| 1566 | C |
---|
| 1567 | C Ecriture etat initial physique |
---|
| 1568 | C |
---|
| 1569 | |
---|
| 1570 | call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm, |
---|
| 1571 | & nqtot,dtphys,real(day_ini),0.0, |
---|
| 1572 | & airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe) |
---|
| 1573 | call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nq, |
---|
| 1574 | & dtphys,real(day_ini), |
---|
| 1575 | & tsurf,tsoil,emis,q2,qsurf, |
---|
| 1576 | & cloudfrac,totalfrac,hice, |
---|
| 1577 | & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) |
---|
| 1578 | |
---|
| 1579 | c======================================================================= |
---|
| 1580 | c Formats |
---|
| 1581 | c======================================================================= |
---|
| 1582 | |
---|
| 1583 | 1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dema |
---|
| 1584 | *rrage est differente de la valeur parametree iim =',i4//) |
---|
| 1585 | 2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dema |
---|
| 1586 | *rrage est differente de la valeur parametree jjm =',i4//) |
---|
| 1587 | 3 FORMAT(//10x,'la valeur de lllm =',i4,2x,'lue sur le fichier demar |
---|
| 1588 | *rage est differente de la valeur parametree llm =',i4//) |
---|
| 1589 | |
---|
| 1590 | write(*,*) "newstart: All is well that ends well." |
---|
| 1591 | |
---|
| 1592 | end |
---|
| 1593 | |
---|