Changeset 1524 for trunk/LMDZ.GENERIC/libf
- Timestamp:
- Mar 29, 2016, 11:45:49 AM (9 years ago)
- Location:
- trunk/LMDZ.GENERIC/libf
- Files:
-
- 1 added
- 1 deleted
- 15 edited
- 1 moved
-
dynlonlat_phylonlat/phystd/iniphysiq_mod.F90 (modified) (1 diff)
-
dynlonlat_phylonlat/phystd/newstart.F (modified) (2 diffs)
-
phystd/comcstfi_mod.F90 (modified) (1 diff)
-
phystd/dyn1d/rcm1d.F (modified) (4 diffs)
-
phystd/dyn1d/temps_mod.F90 (deleted)
-
phystd/inifis_mod.F90 (moved) (moved from trunk/LMDZ.GENERIC/libf/phystd/inifis.F) (6 diffs)
-
phystd/iniwrite.F (modified) (1 diff)
-
phystd/iniwrite_specIR.F (modified) (1 diff)
-
phystd/iniwrite_specVI.F (modified) (1 diff)
-
phystd/phyredem.F90 (modified) (1 diff)
-
phystd/physiq.F90 (modified) (42 diffs)
-
phystd/soil.F (modified) (1 diff)
-
phystd/tabfi.F (modified) (2 diffs)
-
phystd/time_phylmdz_mod.F90 (added)
-
phystd/writediagfi.F (modified) (1 diff)
-
phystd/writediagspecIR.F (modified) (1 diff)
-
phystd/writediagspecVI.F (modified) (1 diff)
-
phystd/writeg1d.F (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90
r1523 r1524 24 24 use planete_mod, only: ini_planete_mod 25 25 USE comvert_mod, ONLY: ap,bp,preff 26 use inifis_mod, only: inifis 26 27 use regular_lonlat_mod, only: init_regular_lonlat, & 27 28 east, west, north, south, & -
trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/newstart.F
r1470 r1524 22 22 USE comgeomfi_h, ONLY: lati, long, area 23 23 use datafile_mod, only: datadir, surfdir 24 ! to use 'getin' 25 ! USE ioipsl_getincom, only: getin 26 USE ioipsl_getincom_p, only: getin_p 24 use ioipsl_getin_p_mod, only: getin_p 27 25 use control_mod, only: day_step, iphysiq, anneeref, planet_type 28 26 use phyredem, only: physdem0, physdem1 … … 37 35 USE temps_mod, ONLY: day_ini 38 36 USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 37 use inifis_mod, only: inifis 39 38 implicit none 40 39 -
trunk/LMDZ.GENERIC/libf/phystd/comcstfi_mod.F90
r1520 r1524 8 8 REAL,SAVE :: cpp ! Cp of the atmosphere 9 9 REAL,SAVE :: rcp ! r/cpp 10 REAL,SAVE :: dtphys ! physics time step (s)11 REAL,SAVE :: daysec ! length of day (s)12 10 REAL,SAVE :: mugaz ! molar mass of the atmosphere (g/mol) 13 11 REAL,SAVE :: omeg ! planet rotation rate (rad/s) 14 12 REAL,SAVE :: avocado ! something like 6.022e23 15 !$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp, dtphys,daysec,mugaz,omeg,avocado)13 !$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,mugaz,omeg,avocado) 16 14 17 15 END MODULE comcstfi_mod -
trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F
r1494 r1524 19 19 & obliquit,nres,z0,lmixmin,emin_turb,coefvis,coefir, 20 20 & timeperi,e_elips,p_elips 21 use comcstfi_mod, only: pi, cpp, daysec, dtphys,rad, g, r,21 use comcstfi_mod, only: pi, cpp, rad, g, r, 22 22 & mugaz, rcp, omeg 23 use time_phylmdz_mod, only: daysec, dtphys 23 24 use callkeys_mod, only: tracer,check_cpp_match,rings_shadow, 24 & 25 & specOLR,water,pceil,ok_slab_ocean 25 26 USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff 26 27 USE logic_mod, ONLY: hybrid,autozlevs 28 use inifis_mod, only: inifis 27 29 implicit none 28 30 … … 130 132 character*20,allocatable :: nametrac(:) ! name of the tracer (no need for adv trac common) 131 133 132 real :: latitude , longitude134 real :: latitude(1), longitude(1) 133 135 134 136 c======================================================================= … … 339 341 longitude=longitude*pi/180.E+0 340 342 343 344 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 345 !!!! PLANETARY CONSTANTS !!!! 346 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 347 348 g = -99999. 349 PRINT *,'GRAVITY in m s-2 ?' 350 call getin("g",g) 351 IF (g.eq.-99999.) THEN 352 PRINT *,"STOP. I NEED g IN RCM1D.DEF." 353 STOP 354 ELSE 355 PRINT *,"--> g = ",g 356 ENDIF 357 358 rad = -99999. 359 PRINT *,'PLANETARY RADIUS in m ?' 360 call getin("rad",rad) 361 ! Planetary radius is needed to compute shadow of the rings 362 IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN 363 PRINT *,"STOP. I NEED rad IN RCM1D.DEF." 364 STOP 365 ELSE 366 PRINT *,"--> rad = ",rad 367 ENDIF 368 369 daysec = -99999. 370 PRINT *,'LENGTH OF A DAY in s ?' 371 call getin("daysec",daysec) 372 IF (daysec.eq.-99999.) THEN 373 PRINT *,"STOP. I NEED daysec IN RCM1D.DEF." 374 STOP 375 ELSE 376 PRINT *,"--> daysec = ",daysec 377 ENDIF 378 omeg=4.*asin(1.)/(daysec) 379 PRINT *,"OK. FROM THIS I WORKED OUT:" 380 PRINT *,"--> omeg = ",omeg 381 382 year_day = -99999. 383 PRINT *,'LENGTH OF A YEAR in days ?' 384 call getin("year_day",year_day) 385 IF (year_day.eq.-99999.) THEN 386 PRINT *,"STOP. I NEED year_day IN RCM1D.DEF." 387 STOP 388 ELSE 389 PRINT *,"--> year_day = ",year_day 390 ENDIF 391 392 periastr = -99999. 393 PRINT *,'MIN DIST STAR-PLANET in AU ?' 394 call getin("periastr",periastr) 395 IF (periastr.eq.-99999.) THEN 396 PRINT *,"STOP. I NEED periastr IN RCM1D.DEF." 397 STOP 398 ELSE 399 PRINT *,"--> periastr = ",periastr 400 ENDIF 401 402 apoastr = -99999. 403 PRINT *,'MAX DIST STAR-PLANET in AU ?' 404 call getin("apoastr",apoastr) 405 IF (apoastr.eq.-99999.) THEN 406 PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF." 407 STOP 408 ELSE 409 PRINT *,"--> apoastr = ",apoastr 410 ENDIF 411 412 peri_day = -99999. 413 PRINT *,'DATE OF PERIASTRON in days ?' 414 call getin("peri_day",peri_day) 415 IF (peri_day.eq.-99999.) THEN 416 PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF." 417 STOP 418 ELSE IF (peri_day.gt.year_day) THEN 419 PRINT *,"STOP. peri_day.gt.year_day" 420 STOP 421 ELSE 422 PRINT *,"--> peri_day = ", peri_day 423 ENDIF 424 425 obliquit = -99999. 426 PRINT *,'OBLIQUITY in deg ?' 427 call getin("obliquit",obliquit) 428 IF (obliquit.eq.-99999.) THEN 429 PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF." 430 STOP 431 ELSE 432 PRINT *,"--> obliquit = ",obliquit 433 ENDIF 434 435 psurf = -99999. 436 PRINT *,'SURFACE PRESSURE in Pa ?' 437 call getin("psurf",psurf) 438 IF (psurf.eq.-99999.) THEN 439 PRINT *,"STOP. I NEED psurf IN RCM1D.DEF." 440 STOP 441 ELSE 442 PRINT *,"--> psurf = ",psurf 443 ENDIF 444 !! we need reference pressures 445 pa=psurf/30. 446 preff=psurf ! these values are not needed in 1D (are you sure JL12) 447 448 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 449 !!!! END PLANETARY CONSTANTS !!!! 450 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 451 452 c Date et heure locale du debut du run 453 c ------------------------------------ 454 c Date (en sols depuis le solstice de printemps) du debut du run 455 day0 = 0 ! default value for day0 456 write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?' 457 call getin("day0",day0) 458 day=float(day0) 459 write(*,*) " day0 = ",day0 460 c Heure de demarrage 461 time=0 ! default value for time 462 write(*,*)'Initial local time (in hours, between 0 and 24)?' 463 call getin("time",time) 464 write(*,*)" time = ",time 465 time=time/24.E+0 ! convert time (hours) to fraction of sol 466 467 468 c Discretisation (Definition de la grille et des pas de temps) 469 c -------------- 470 c 471 nlayer=llm 472 nlevel=nlayer+1 473 nsoil=nsoilmx 474 475 day_step=48 ! default value for day_step 476 PRINT *,'Number of time steps per sol ?' 477 call getin("day_step",day_step) 478 write(*,*) " day_step = ",day_step 479 480 481 ecritphy=day_step ! default value for ecritphy 482 PRINT *,'Nunber of steps between writediagfi ?' 483 call getin("ecritphy",ecritphy) 484 write(*,*) " ecritphy = ",ecritphy 485 486 ndt=10 ! default value for ndt 487 PRINT *,'Number of sols to run ?' 488 call getin("ndt",ndt) 489 write(*,*) " ndt = ",ndt 490 491 ndt=ndt*day_step 492 dtphys=daysec/day_step 493 write(*,*)"-------------------------------------" 494 write(*,*)"-------------------------------------" 495 write(*,*)"--> Physical timestep is ",dtphys 496 write(*,*)"-------------------------------------" 497 write(*,*)"-------------------------------------" 498 341 499 !!! CALL INIFIS 342 500 !!! - read callphys.def … … 355 513 r = 8.314511E+0 * 1000.E+0 / mugaz 356 514 rcp = r / cpp 357 358 359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!360 !!!! PLANETARY CONSTANTS !!!!361 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!362 363 g = -99999.364 PRINT *,'GRAVITY in m s-2 ?'365 call getin("g",g)366 IF (g.eq.-99999.) THEN367 PRINT *,"STOP. I NEED g IN RCM1D.DEF."368 STOP369 ELSE370 PRINT *,"--> g = ",g371 ENDIF372 373 rad = -99999.374 PRINT *,'PLANETARY RADIUS in m ?'375 call getin("rad",rad)376 ! Planetary radius is needed to compute shadow of the rings377 IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN378 PRINT *,"STOP. I NEED rad IN RCM1D.DEF."379 STOP380 ELSE381 PRINT *,"--> rad = ",rad382 ENDIF383 384 daysec = -99999.385 PRINT *,'LENGTH OF A DAY in s ?'386 call getin("daysec",daysec)387 IF (daysec.eq.-99999.) THEN388 PRINT *,"STOP. I NEED daysec IN RCM1D.DEF."389 STOP390 ELSE391 PRINT *,"--> daysec = ",daysec392 ENDIF393 omeg=4.*asin(1.)/(daysec)394 PRINT *,"OK. FROM THIS I WORKED OUT:"395 PRINT *,"--> omeg = ",omeg396 397 year_day = -99999.398 PRINT *,'LENGTH OF A YEAR in days ?'399 call getin("year_day",year_day)400 IF (year_day.eq.-99999.) THEN401 PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."402 STOP403 ELSE404 PRINT *,"--> year_day = ",year_day405 ENDIF406 407 periastr = -99999.408 PRINT *,'MIN DIST STAR-PLANET in AU ?'409 call getin("periastr",periastr)410 IF (periastr.eq.-99999.) THEN411 PRINT *,"STOP. I NEED periastr IN RCM1D.DEF."412 STOP413 ELSE414 PRINT *,"--> periastr = ",periastr415 ENDIF416 417 apoastr = -99999.418 PRINT *,'MAX DIST STAR-PLANET in AU ?'419 call getin("apoastr",apoastr)420 IF (apoastr.eq.-99999.) THEN421 PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF."422 STOP423 ELSE424 PRINT *,"--> apoastr = ",apoastr425 ENDIF426 427 peri_day = -99999.428 PRINT *,'DATE OF PERIASTRON in days ?'429 call getin("peri_day",peri_day)430 IF (peri_day.eq.-99999.) THEN431 PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF."432 STOP433 ELSE IF (peri_day.gt.year_day) THEN434 PRINT *,"STOP. peri_day.gt.year_day"435 STOP436 ELSE437 PRINT *,"--> peri_day = ", peri_day438 ENDIF439 440 obliquit = -99999.441 PRINT *,'OBLIQUITY in deg ?'442 call getin("obliquit",obliquit)443 IF (obliquit.eq.-99999.) THEN444 PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF."445 STOP446 ELSE447 PRINT *,"--> obliquit = ",obliquit448 ENDIF449 450 psurf = -99999.451 PRINT *,'SURFACE PRESSURE in Pa ?'452 call getin("psurf",psurf)453 IF (psurf.eq.-99999.) THEN454 PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."455 STOP456 ELSE457 PRINT *,"--> psurf = ",psurf458 ENDIF459 !! we need reference pressures460 pa=psurf/30.461 preff=psurf ! these values are not needed in 1D (are you sure JL12)462 463 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!464 !!!! END PLANETARY CONSTANTS !!!!465 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!466 467 c Date et heure locale du debut du run468 c ------------------------------------469 c Date (en sols depuis le solstice de printemps) du debut du run470 day0 = 0 ! default value for day0471 write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'472 call getin("day0",day0)473 day=float(day0)474 write(*,*) " day0 = ",day0475 c Heure de demarrage476 time=0 ! default value for time477 write(*,*)'Initial local time (in hours, between 0 and 24)?'478 call getin("time",time)479 write(*,*)" time = ",time480 time=time/24.E+0 ! convert time (hours) to fraction of sol481 482 483 c Discretisation (Definition de la grille et des pas de temps)484 c --------------485 c486 nlayer=llm487 nlevel=nlayer+1488 nsoil=nsoilmx489 490 day_step=48 ! default value for day_step491 PRINT *,'Number of time steps per sol ?'492 call getin("day_step",day_step)493 write(*,*) " day_step = ",day_step494 495 496 ecritphy=day_step ! default value for ecritphy497 PRINT *,'Nunber of steps between writediagfi ?'498 call getin("ecritphy",ecritphy)499 write(*,*) " ecritphy = ",ecritphy500 501 ndt=10 ! default value for ndt502 PRINT *,'Number of sols to run ?'503 call getin("ndt",ndt)504 write(*,*) " ndt = ",ndt505 506 ndt=ndt*day_step507 dtphys=daysec/day_step508 write(*,*)"-------------------------------------"509 write(*,*)"-------------------------------------"510 write(*,*)"--> Physical timestep is ",dtphys511 write(*,*)"-------------------------------------"512 write(*,*)"-------------------------------------"513 515 514 516 c output spectrum? -
trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
r1521 r1524 1 SUBROUTINE inifis(ngrid,nlayer,nq, 2 $ day_ini,pdaysec,ptimestep, 3 $ plat,plon,parea, 4 $ prad,pg,pr,pcpp) 5 6 use radinc_h, only : naerkind 7 use datafile_mod, only: datadir 8 use comdiurn_h, only: sinlat, coslat, sinlon, coslon 9 use comgeomfi_h, only: long, lati, area, totarea, totarea_planet 10 use comsoil_h, only: ini_comsoil_h 11 use control_mod, only: ecritphy 12 use comcstfi_mod, only: rad, daysec, dtphys, cpp, g, r, rcp, 13 & mugaz, pi, avocado 14 use planete_mod, only: nres 15 use planetwide_mod, only: planetwide_sumval 16 use callkeys_mod 1 MODULE inifis_mod 2 IMPLICIT NONE 3 4 CONTAINS 5 6 SUBROUTINE inifis(ngrid,nlayer,nq, & 7 day_ini,pdaysec,ptimestep, & 8 plat,plon,parea, & 9 prad,pg,pr,pcpp) 10 11 use radinc_h, only : naerkind 12 use datafile_mod, only: datadir 13 use comdiurn_h, only: sinlat, coslat, sinlon, coslon 14 use comgeomfi_h, only: long, lati, area, totarea, totarea_planet 15 use comsoil_h, only: ini_comsoil_h 16 use control_mod, only: ecritphy 17 use comcstfi_mod, only: rad, cpp, g, r, rcp, & 18 mugaz, pi, avocado 19 use planete_mod, only: nres 20 use planetwide_mod, only: planetwide_sumval 21 use callkeys_mod 22 use time_phylmdz_mod, only: init_time, daysec, dtphys 17 23 18 24 !======================================================================= … … 51 57 ! declarations: 52 58 ! ------------- 53 54 55 56 57 58 59 59 use datafile_mod, only: datadir 60 use ioipsl_getin_p_mod, only: getin_p 61 IMPLICIT NONE 62 63 64 65 REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep 60 66 61 62 63 integerday_ini64 67 INTEGER,INTENT(IN) :: ngrid,nlayer,nq 68 REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid) 69 integer,intent(in) :: day_ini 70 INTEGER ig,ierr 65 71 66 67 68 72 EXTERNAL iniorbit,orbite 73 EXTERNAL SSUM 74 REAL SSUM 69 75 70 71 72 73 74 logical :: parameter, doubleq=.false. 75 76 real psurf,pN2 ! added by RW for Gliese 581d N2+CO2 77 78 79 daysec=pdaysec80 dtphys=ptimestep81 cpp=pcpp82 g=pg83 r=pr84 rcp=r/cpp85 86 avocado = 6.02214179e23 ! added by RW87 88 89 90 91 76 CHARACTER ch1*12 77 CHARACTER ch80*80 78 79 logical chem, h2o 80 81 real psurf,pN2 ! added by RW for Gliese 581d N2+CO2 82 83 ! initialize constants in comcstfi_mod 84 rad=prad 85 cpp=pcpp 86 g=pg 87 r=pr 88 rcp=r/cpp 89 pi=2.*asin(1.) 90 avocado = 6.02214179e23 ! added by RW 91 92 ! Initialize some "temporal and calendar" related variables 93 CALL init_time(day_ini,pdaysec,ptimestep) 94 95 ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps) 96 ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON) 97 call getin_p("ecritphy",ecritphy) 92 98 93 99 ! -------------------------------------------------------------- … … 95 101 ! -------------------------------------------------------------- 96 102 97 ! check that 'callphys.def' file is around 98 OPEN(99,file='callphys.def',status='old',form='formatted' 99 & ,iostat=ierr) 100 CLOSE(99) 101 IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module 103 ! check that 'callphys.def' file is around 104 OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr) 105 CLOSE(99) 106 IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module 102 107 103 108 !!! IF(ierr.EQ.0) THEN 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 write(*,*) "Run with or without atm mass update ",122 &" due to tracer evaporation/condensation?"123 124 125 126 127 128 129 130 131 132 133 134 write(*,*) "(if season=false, Ls stays constant, to value ",135 &"set in 'start'"136 137 138 139 140 141 142 143 144 145 146 147 148 109 IF(iscallphys) THEN 110 PRINT* 111 PRINT* 112 PRINT*,'--------------------------------------------' 113 PRINT*,' inifis: Parametres pour la physique (callphys.def)' 114 PRINT*,'--------------------------------------------' 115 116 write(*,*) "Directory where external input files are:" 117 ! default 'datadir' is set in "datadir_mod" 118 call getin_p("datadir",datadir) ! default path 119 write(*,*) " datadir = ",trim(datadir) 120 121 write(*,*) "Run with or without tracer transport ?" 122 tracer=.false. ! default value 123 call getin_p("tracer",tracer) 124 write(*,*) " tracer = ",tracer 125 126 write(*,*) "Run with or without atm mass update ", & 127 " due to tracer evaporation/condensation?" 128 mass_redistrib=.false. ! default value 129 call getin_p("mass_redistrib",mass_redistrib) 130 write(*,*) " mass_redistrib = ",mass_redistrib 131 132 write(*,*) "Diurnal cycle ?" 133 write(*,*) "(if diurnal=false, diurnal averaged solar heating)" 134 diurnal=.true. ! default value 135 call getin_p("diurnal",diurnal) 136 write(*,*) " diurnal = ",diurnal 137 138 write(*,*) "Seasonal cycle ?" 139 write(*,*) "(if season=false, Ls stays constant, to value ", & 140 "set in 'start'" 141 season=.true. ! default value 142 call getin_p("season",season) 143 write(*,*) " season = ",season 144 145 write(*,*) "Tidally resonant rotation ?" 146 tlocked=.false. ! default value 147 call getin_p("tlocked",tlocked) 148 write(*,*) "tlocked = ",tlocked 149 150 write(*,*) "Saturn ring shadowing ?" 151 rings_shadow = .false. 152 call getin_p("rings_shadow", rings_shadow) 153 write(*,*) "rings_shadow = ", rings_shadow 149 154 150 151 152 153 154 155 156 157 158 155 write(*,*) "Compute latitude-dependent gravity field?" 156 oblate = .false. 157 call getin_p("oblate", oblate) 158 write(*,*) "oblate = ", oblate 159 160 write(*,*) "Flattening of the planet (a-b)/a " 161 flatten = 0.0 162 call getin_p("flatten", flatten) 163 write(*,*) "flatten = ", flatten 159 164 160 165 161 162 163 164 166 write(*,*) "Needed if oblate=.true.: J2" 167 J2 = 0.0 168 call getin_p("J2", J2) 169 write(*,*) "J2 = ", J2 165 170 166 167 168 169 170 171 172 173 174 171 write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)" 172 MassPlanet = 0.0 173 call getin_p("MassPlanet", MassPlanet) 174 write(*,*) "MassPlanet = ", MassPlanet 175 176 write(*,*) "Needed if oblate=.true.: Planet mean radius (m)" 177 Rmean = 0.0 178 call getin_p("Rmean", Rmean) 179 write(*,*) "Rmean = ", Rmean 175 180 176 181 ! Test of incompatibility: 177 182 ! if tlocked, then diurnal should be false 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 write(*,*) "call gaseous absorption in the visible bands?",224 &"(matters only if callrad=T)"225 226 227 183 if (tlocked.and.diurnal) then 184 print*,'If diurnal=true, we should turn off tlocked.' 185 stop 186 endif 187 188 write(*,*) "Tidal resonance ratio ?" 189 nres=0 ! default value 190 call getin_p("nres",nres) 191 write(*,*) "nres = ",nres 192 193 write(*,*) "Write some extra output to the screen ?" 194 lwrite=.false. ! default value 195 call getin_p("lwrite",lwrite) 196 write(*,*) " lwrite = ",lwrite 197 198 write(*,*) "Save statistics in file stats.nc ?" 199 callstats=.true. ! default value 200 call getin_p("callstats",callstats) 201 write(*,*) " callstats = ",callstats 202 203 write(*,*) "Test energy conservation of model physics ?" 204 enertest=.false. ! default value 205 call getin_p("enertest",enertest) 206 write(*,*) " enertest = ",enertest 207 208 write(*,*) "Check to see if cpp values used match gases.def ?" 209 check_cpp_match=.true. ! default value 210 call getin_p("check_cpp_match",check_cpp_match) 211 write(*,*) " check_cpp_match = ",check_cpp_match 212 213 write(*,*) "call radiative transfer ?" 214 callrad=.true. ! default value 215 call getin_p("callrad",callrad) 216 write(*,*) " callrad = ",callrad 217 218 write(*,*) "call correlated-k radiative transfer ?" 219 corrk=.true. ! default value 220 call getin_p("corrk",corrk) 221 write(*,*) " corrk = ",corrk 222 223 write(*,*) "prohibit calculations outside corrk T grid?" 224 strictboundcorrk=.true. ! default value 225 call getin_p("strictboundcorrk",strictboundcorrk) 226 write(*,*) "strictboundcorrk = ",strictboundcorrk 227 228 write(*,*) "call gaseous absorption in the visible bands?", & 229 "(matters only if callrad=T)" 230 callgasvis=.false. ! default value 231 call getin_p("callgasvis",callgasvis) 232 write(*,*) " callgasvis = ",callgasvis 228 233 229 write(*,*) "call continuum opacities in radiative transfer ?",230 &"(matters only if callrad=T)"231 232 233 234 235 236 237 238 234 write(*,*) "call continuum opacities in radiative transfer ?", & 235 "(matters only if callrad=T)" 236 continuum=.true. ! default value 237 call getin_p("continuum",continuum) 238 write(*,*) " continuum = ",continuum 239 240 write(*,*) "use analytic function for H2O continuum ?" 241 H2Ocont_simple=.false. ! default value 242 call getin_p("H2Ocont_simple",H2Ocont_simple) 243 write(*,*) " H2Ocont_simple = ",H2Ocont_simple 239 244 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 245 write(*,*) "call turbulent vertical diffusion ?" 246 calldifv=.true. ! default value 247 call getin_p("calldifv",calldifv) 248 write(*,*) " calldifv = ",calldifv 249 250 write(*,*) "use turbdiff instead of vdifc ?" 251 UseTurbDiff=.true. ! default value 252 call getin_p("UseTurbDiff",UseTurbDiff) 253 write(*,*) " UseTurbDiff = ",UseTurbDiff 254 255 write(*,*) "call convective adjustment ?" 256 calladj=.true. ! default value 257 call getin_p("calladj",calladj) 258 write(*,*) " calladj = ",calladj 259 260 write(*,*) "call CO2 condensation ?" 261 co2cond=.false. ! default value 262 call getin_p("co2cond",co2cond) 263 write(*,*) " co2cond = ",co2cond 259 264 ! Test of incompatibility 260 261 262 263 265 if (co2cond.and.(.not.tracer)) then 266 print*,'We need a CO2 ice tracer to condense CO2' 267 call abort 268 endif 264 269 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 270 write(*,*) "CO2 supersaturation level ?" 271 co2supsat=1.0 ! default value 272 call getin_p("co2supsat",co2supsat) 273 write(*,*) " co2supsat = ",co2supsat 274 275 write(*,*) "Radiative timescale for Newtonian cooling ?" 276 tau_relax=30. ! default value 277 call getin_p("tau_relax",tau_relax) 278 write(*,*) " tau_relax = ",tau_relax 279 tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds 280 281 write(*,*)"call thermal conduction in the soil ?" 282 callsoil=.true. ! default value 283 call getin_p("callsoil",callsoil) 284 write(*,*) " callsoil = ",callsoil 280 285 281 write(*,*)"Rad transfer is computed every iradia",282 &" physical timestep"283 284 285 286 write(*,*)"Rad transfer is computed every iradia", & 287 " physical timestep" 288 iradia=1 ! default value 289 call getin_p("iradia",iradia) 290 write(*,*)" iradia = ",iradia 286 291 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 292 write(*,*)"Rayleigh scattering ?" 293 rayleigh=.false. 294 call getin_p("rayleigh",rayleigh) 295 write(*,*)" rayleigh = ",rayleigh 296 297 write(*,*) "Use blackbody for stellar spectrum ?" 298 stelbbody=.false. ! default value 299 call getin_p("stelbbody",stelbbody) 300 write(*,*) " stelbbody = ",stelbbody 301 302 write(*,*) "Stellar blackbody temperature ?" 303 stelTbb=5800.0 ! default value 304 call getin_p("stelTbb",stelTbb) 305 write(*,*) " stelTbb = ",stelTbb 306 307 write(*,*)"Output mean OLR in 1D?" 308 meanOLR=.false. 309 call getin_p("meanOLR",meanOLR) 310 write(*,*)" meanOLR = ",meanOLR 311 312 write(*,*)"Output spectral OLR in 3D?" 313 specOLR=.false. 314 call getin_p("specOLR",specOLR) 315 write(*,*)" specOLR = ",specOLR 316 317 write(*,*)"Operate in kastprof mode?" 318 kastprof=.false. 319 call getin_p("kastprof",kastprof) 320 write(*,*)" kastprof = ",kastprof 321 322 write(*,*)"Uniform absorption in radiative transfer?" 323 graybody=.false. 324 call getin_p("graybody",graybody) 325 write(*,*)" graybody = ",graybody 321 326 322 327 ! Slab Ocean 323 324 325 326 327 328 329 330 331 332 333 334 335 336 328 write(*,*) "Use slab-ocean ?" 329 ok_slab_ocean=.false. ! default value 330 call getin_p("ok_slab_ocean",ok_slab_ocean) 331 write(*,*) "ok_slab_ocean = ",ok_slab_ocean 332 333 write(*,*) "Use slab-sea-ice ?" 334 ok_slab_sic=.true. ! default value 335 call getin_p("ok_slab_sic",ok_slab_sic) 336 write(*,*) "ok_slab_sic = ",ok_slab_sic 337 338 write(*,*) "Use heat transport for the ocean ?" 339 ok_slab_heat_transp=.true. ! default value 340 call getin_p("ok_slab_heat_transp",ok_slab_heat_transp) 341 write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp 337 342 338 343 … … 340 345 ! Test of incompatibility: 341 346 ! if kastprof used, we must be in 1D 342 343 344 345 346 347 348 349 350 351 352 353 354 355 347 if (kastprof.and.(ngrid.gt.1)) then 348 print*,'kastprof can only be used in 1D!' 349 call abort 350 endif 351 352 write(*,*)"Stratospheric temperature for kastprof mode?" 353 Tstrat=167.0 354 call getin_p("Tstrat",Tstrat) 355 write(*,*)" Tstrat = ",Tstrat 356 357 write(*,*)"Remove lower boundary?" 358 nosurf=.false. 359 call getin_p("nosurf",nosurf) 360 write(*,*)" nosurf = ",nosurf 356 361 357 362 ! Tests of incompatibility: 358 359 360 361 362 363 364 write(*,*)"Add an internal heat flux?",365 ."... matters only if callsoil=F"366 367 368 369 370 371 372 373 363 if (nosurf.and.callsoil) then 364 print*,'nosurf not compatible with soil scheme!' 365 print*,'... got to make a choice!' 366 call abort 367 endif 368 369 write(*,*)"Add an internal heat flux?", & 370 "... matters only if callsoil=F" 371 intheat=0. 372 call getin_p("intheat",intheat) 373 write(*,*)" intheat = ",intheat 374 375 write(*,*)"Use Newtonian cooling for radiative transfer?" 376 newtonian=.false. 377 call getin_p("newtonian",newtonian) 378 write(*,*)" newtonian = ",newtonian 374 379 375 380 ! Tests of incompatibility: 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 381 if (newtonian.and.corrk) then 382 print*,'newtonian not compatible with correlated-k!' 383 call abort 384 endif 385 if (newtonian.and.calladj) then 386 print*,'newtonian not compatible with adjustment!' 387 call abort 388 endif 389 if (newtonian.and.calldifv) then 390 print*,'newtonian not compatible with a boundary layer!' 391 call abort 392 endif 393 394 write(*,*)"Test physics timescale in 1D?" 395 testradtimes=.false. 396 call getin_p("testradtimes",testradtimes) 397 write(*,*)" testradtimes = ",testradtimes 393 398 394 399 ! Test of incompatibility: 395 400 ! if testradtimes used, we must be in 1D 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 401 if (testradtimes.and.(ngrid.gt.1)) then 402 print*,'testradtimes can only be used in 1D!' 403 call abort 404 endif 405 406 write(*,*)"Default planetary temperature?" 407 tplanet=215.0 408 call getin_p("tplanet",tplanet) 409 write(*,*)" tplanet = ",tplanet 410 411 write(*,*)"Which star?" 412 startype=1 ! default value = Sol 413 call getin_p("startype",startype) 414 write(*,*)" startype = ",startype 415 416 write(*,*)"Value of stellar flux at 1 AU?" 417 Fat1AU=1356.0 ! default value = Sol today 418 call getin_p("Fat1AU",Fat1AU) 419 write(*,*)" Fat1AU = ",Fat1AU 415 420 416 421 417 422 ! TRACERS: 418 423 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 424 write(*,*)"Varying H2O cloud fraction?" 425 CLFvarying=.false. ! default value 426 call getin_p("CLFvarying",CLFvarying) 427 write(*,*)" CLFvarying = ",CLFvarying 428 429 write(*,*)"Value of fixed H2O cloud fraction?" 430 CLFfixval=1.0 ! default value 431 call getin_p("CLFfixval",CLFfixval) 432 write(*,*)" CLFfixval = ",CLFfixval 433 434 write(*,*)"fixed radii for Cloud particles?" 435 radfixed=.false. ! default value 436 call getin_p("radfixed",radfixed) 437 write(*,*)" radfixed = ",radfixed 438 439 if(kastprof)then 440 radfixed=.true. 441 endif 442 443 write(*,*)"Number mixing ratio of CO2 ice particles:" 444 Nmix_co2=1.e6 ! default value 445 call getin_p("Nmix_co2",Nmix_co2) 446 write(*,*)" Nmix_co2 = ",Nmix_co2 442 447 443 448 ! write(*,*)"Number of radiatively active aerosols:" … … 446 451 ! write(*,*)" naerkind = ",naerkind 447 452 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 453 write(*,*)"Opacity of dust (if used):" 454 dusttau=0. ! default value 455 call getin_p("dusttau",dusttau) 456 write(*,*)" dusttau = ",dusttau 457 458 write(*,*)"Radiatively active CO2 aerosols?" 459 aeroco2=.false. ! default value 460 call getin_p("aeroco2",aeroco2) 461 write(*,*)" aeroco2 = ",aeroco2 462 463 write(*,*)"Fixed CO2 aerosol distribution?" 464 aerofixco2=.false. ! default value 465 call getin_p("aerofixco2",aerofixco2) 466 write(*,*)" aerofixco2 = ",aerofixco2 467 468 write(*,*)"Radiatively active water ice?" 469 aeroh2o=.false. ! default value 470 call getin_p("aeroh2o",aeroh2o) 471 write(*,*)" aeroh2o = ",aeroh2o 472 473 write(*,*)"Fixed H2O aerosol distribution?" 474 aerofixh2o=.false. ! default value 475 call getin_p("aerofixh2o",aerofixh2o) 476 write(*,*)" aerofixh2o = ",aerofixh2o 477 478 write(*,*)"Radiatively active sulfuric acid aersols?" 479 aeroh2so4=.false. ! default value 480 call getin_p("aeroh2so4",aeroh2so4) 481 write(*,*)" aeroh2so4 = ",aeroh2so4 477 482 478 483 !================================= 479 484 480 481 482 483 484 485 write(*,*)"TWOLAY AEROSOL: total optical depth ",486 &"in the tropospheric layer (visible)"487 488 489 490 491 write(*,*)"TWOLAY AEROSOL: total optical depth ",492 &"in the stratospheric layer (visible)"493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 write(*,*)"TWOLAY AEROSOL: particle size in the ",518 &"tropospheric layer, in meters"519 520 521 522 523 write(*,*)"TWOLAY AEROSOL: particle size in the ",524 &"stratospheric layer, in meters"525 526 527 485 write(*,*)"Radiatively active two-layer aersols?" 486 aeroback2lay=.false. ! default value 487 call getin_p("aeroback2lay",aeroback2lay) 488 write(*,*)" aeroback2lay = ",aeroback2lay 489 490 write(*,*)"TWOLAY AEROSOL: total optical depth ", & 491 "in the tropospheric layer (visible)" 492 obs_tau_col_tropo=8.D0 493 call getin_p("obs_tau_col_tropo",obs_tau_col_tropo) 494 write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo 495 496 write(*,*)"TWOLAY AEROSOL: total optical depth ", & 497 "in the stratospheric layer (visible)" 498 obs_tau_col_strato=0.08D0 499 call getin_p("obs_tau_col_strato",obs_tau_col_strato) 500 write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato 501 502 write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa" 503 pres_bottom_tropo=66000.0 504 call getin_p("pres_bottom_tropo",pres_bottom_tropo) 505 write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo 506 507 write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa" 508 pres_top_tropo=18000.0 509 call getin_p("pres_top_tropo",pres_top_tropo) 510 write(*,*)" pres_top_tropo = ",pres_top_tropo 511 512 write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa" 513 pres_bottom_strato=2000.0 514 call getin_p("pres_bottom_strato",pres_bottom_strato) 515 write(*,*)" pres_bottom_strato = ",pres_bottom_strato 516 517 write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa" 518 pres_top_strato=100.0 519 call getin_p("pres_top_strato",pres_top_strato) 520 write(*,*)" pres_top_strato = ",pres_top_strato 521 522 write(*,*)"TWOLAY AEROSOL: particle size in the ", & 523 "tropospheric layer, in meters" 524 size_tropo=2.e-6 525 call getin_p("size_tropo",size_tropo) 526 write(*,*)" size_tropo = ",size_tropo 527 528 write(*,*)"TWOLAY AEROSOL: particle size in the ", & 529 "stratospheric layer, in meters" 530 size_strato=1.e-7 531 call getin_p("size_strato",size_strato) 532 write(*,*)" size_strato = ",size_strato 528 533 529 534 !================================= 530 535 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 536 write(*,*)"Cloud pressure level (with kastprof only):" 537 cloudlvl=0. ! default value 538 call getin_p("cloudlvl",cloudlvl) 539 write(*,*)" cloudlvl = ",cloudlvl 540 541 write(*,*)"Is the variable gas species radiatively active?" 542 Tstrat=167.0 543 varactive=.false. 544 call getin_p("varactive",varactive) 545 write(*,*)" varactive = ",varactive 546 547 write(*,*)"Is the variable gas species distribution set?" 548 varfixed=.false. 549 call getin_p("varfixed",varfixed) 550 write(*,*)" varfixed = ",varfixed 551 552 write(*,*)"What is the saturation % of the variable species?" 553 satval=0.8 554 call getin_p("satval",satval) 555 write(*,*)" satval = ",satval 551 556 552 557 553 558 ! Test of incompatibility: 554 559 ! if varactive, then varfixed should be false 555 556 557 558 559 560 561 562 563 564 565 566 567 568 560 if (varactive.and.varfixed) then 561 print*,'if varactive, varfixed must be OFF!' 562 stop 563 endif 564 565 write(*,*) "Gravitationnal sedimentation ?" 566 sedimentation=.false. ! default value 567 call getin_p("sedimentation",sedimentation) 568 write(*,*) " sedimentation = ",sedimentation 569 570 write(*,*) "Compute water cycle ?" 571 water=.false. ! default value 572 call getin_p("water",water) 573 write(*,*) " water = ",water 569 574 570 575 ! Test of incompatibility: 571 576 ! if water is true, there should be at least a tracer 572 573 574 575 576 577 578 579 580 577 if (water.and.(.not.tracer)) then 578 print*,'if water is ON, tracer must be ON too!' 579 stop 580 endif 581 582 write(*,*) "Include water condensation ?" 583 watercond=.false. ! default value 584 call getin_p("watercond",watercond) 585 write(*,*) " watercond = ",watercond 581 586 582 587 ! Test of incompatibility: 583 588 ! if watercond is used, then water should be used too 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 589 if (watercond.and.(.not.water)) then 590 print*,'if watercond is used, water should be used too' 591 stop 592 endif 593 594 write(*,*) "Include water precipitation ?" 595 waterrain=.false. ! default value 596 call getin_p("waterrain",waterrain) 597 write(*,*) " waterrain = ",waterrain 598 599 write(*,*) "Include surface hydrology ?" 600 hydrology=.false. ! default value 601 call getin_p("hydrology",hydrology) 602 write(*,*) " hydrology = ",hydrology 603 604 write(*,*) "Evolve surface water sources ?" 605 sourceevol=.false. ! default value 606 call getin_p("sourceevol",sourceevol) 607 write(*,*) " sourceevol = ",sourceevol 608 609 write(*,*) "Ice evolution timestep ?" 610 icetstep=100.0 ! default value 611 call getin_p("icetstep",icetstep) 612 write(*,*) " icetstep = ",icetstep 608 613 609 610 611 612 613 614 615 616 617 618 619 614 write(*,*) "Spectral Dependant albedo ?" 615 albedo_spectral_mode=.false. ! default value 616 call getin_p("albedo_spectral_mode",albedo_spectral_mode) 617 write(*,*) " albedo_spectral_mode = ",albedo_spectral_mode 618 619 write(*,*) "Snow albedo ?" 620 write(*,*) "If albedo_spectral_mode=.true., then this " 621 write(*,*) "corresponds to the 0.5 microns snow albedo." 622 albedosnow=0.5 ! default value 623 call getin_p("albedosnow",albedosnow) 624 write(*,*) " albedosnow = ",albedosnow 620 625 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 626 write(*,*) "CO2 ice albedo ?" 627 albedoco2ice=0.5 ! default value 628 call getin_p("albedoco2ice",albedoco2ice) 629 write(*,*) " albedoco2ice = ",albedoco2ice 630 631 write(*,*) "Maximum ice thickness ?" 632 maxicethick=2.0 ! default value 633 call getin_p("maxicethick",maxicethick) 634 write(*,*) " maxicethick = ",maxicethick 635 636 write(*,*) "Freezing point of seawater ?" 637 Tsaldiff=-1.8 ! default value 638 call getin_p("Tsaldiff",Tsaldiff) 639 write(*,*) " Tsaldiff = ",Tsaldiff 640 641 write(*,*) "Does user want to force cpp and mugaz?" 642 force_cpp=.false. ! default value 643 call getin_p("force_cpp",force_cpp) 644 write(*,*) " force_cpp = ",force_cpp 645 646 if (force_cpp) then 647 mugaz = -99999. 648 PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?' 649 call getin_p("mugaz",mugaz) 650 IF (mugaz.eq.-99999.) THEN 651 PRINT *, "mugaz must be set if force_cpp = T" 652 STOP 653 ELSE 654 write(*,*) "mugaz=",mugaz 655 ENDIF 656 cpp = -99999. 657 PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?' 658 call getin_p("cpp",cpp) 659 IF (cpp.eq.-99999.) THEN 660 PRINT *, "cpp must be set if force_cpp = T" 661 STOP 662 ELSE 663 write(*,*) "cpp=",cpp 664 ENDIF 660 665 ! else 661 666 ! mugaz=8.314*1000./pr 662 endif 663 call su_gases 664 call calc_cpp_mugaz 665 666 PRINT*,'--------------------------------------------' 667 PRINT* 668 PRINT* 669 ELSE 670 write(*,*) 671 write(*,*) 'Cannot read file callphys.def. Is it here ?' 672 stop 673 ENDIF 674 675 8000 FORMAT(t5,a12,l8) 676 8001 FORMAT(t5,a12,i8) 677 678 PRINT* 679 PRINT*,'inifis: daysec',daysec 680 PRINT* 681 PRINT*,'inifis: The radiative transfer is computed:' 682 PRINT*,' each ',iradia,' physical time-step' 683 PRINT*,' or each ',iradia*dtphys,' seconds' 684 PRINT* 667 endif ! of if (force_cpp) 668 call su_gases 669 call calc_cpp_mugaz 670 671 PRINT*,'--------------------------------------------' 672 PRINT* 673 PRINT* 674 ELSE 675 write(*,*) 676 write(*,*) 'Cannot read file callphys.def. Is it here ?' 677 stop 678 ENDIF ! of IF(iscallphys) 679 680 PRINT* 681 PRINT*,'inifis: daysec',daysec 682 PRINT* 683 PRINT*,'inifis: The radiative transfer is computed:' 684 PRINT*,' each ',iradia,' physical time-step' 685 PRINT*,' or each ',iradia*dtphys,' seconds' 686 PRINT* 685 687 686 688 … … 689 691 ! ------------------------ 690 692 691 ! ALLOCATE ARRAYS IN comgeomfi_h 692 IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngrid)) 693 IF (.not. ALLOCATED(long)) ALLOCATE(long(ngrid)) 694 IF (.not. ALLOCATED(area)) ALLOCATE(area(ngrid)) 695 696 CALL SCOPY(ngrid,plon,1,long,1) 697 CALL SCOPY(ngrid,plat,1,lati,1) 698 CALL SCOPY(ngrid,parea,1,area,1) 699 totarea=SSUM(ngrid,area,1) 700 call planetwide_sumval(area,totarea_planet) 701 702 !! those are defined in comdiurn_h.F90 703 IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid)) 704 IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid)) 705 IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid)) 706 IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid)) 707 708 DO ig=1,ngrid 709 sinlat(ig)=sin(plat(ig)) 710 coslat(ig)=cos(plat(ig)) 711 sinlon(ig)=sin(plon(ig)) 712 coslon(ig)=cos(plon(ig)) 713 ENDDO 714 715 pi=2.*asin(1.) ! NB: pi is a common in comcstfi_mod 716 717 ! allocate "comsoil_h" arrays 718 call ini_comsoil_h(ngrid) 693 ! ALLOCATE ARRAYS IN comgeomfi_h 694 IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngrid)) 695 IF (.not. ALLOCATED(long)) ALLOCATE(long(ngrid)) 696 IF (.not. ALLOCATED(area)) ALLOCATE(area(ngrid)) 697 698 CALL SCOPY(ngrid,plon,1,long,1) 699 CALL SCOPY(ngrid,plat,1,lati,1) 700 CALL SCOPY(ngrid,parea,1,area,1) 701 totarea=SSUM(ngrid,area,1) 702 call planetwide_sumval(area,totarea_planet) 703 704 !! those are defined in comdiurn_h.F90 705 IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid)) 706 IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid)) 707 IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid)) 708 IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid)) 709 710 DO ig=1,ngrid 711 sinlat(ig)=sin(plat(ig)) 712 coslat(ig)=cos(plat(ig)) 713 sinlon(ig)=sin(plon(ig)) 714 coslon(ig)=cos(plon(ig)) 715 ENDDO 716 717 ! allocate "comsoil_h" arrays 718 call ini_comsoil_h(ngrid) 719 719 720 END 720 END SUBROUTINE inifis 721 722 END MODULE inifis_mod -
trunk/LMDZ.GENERIC/libf/phystd/iniwrite.F
r1422 r1524 2 2 3 3 use comsoil_h, only: mlayer, nsoilmx 4 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, daysec, dtphys,5 & pi4 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi 5 use time_phylmdz_mod, only: daysec, dtphys 6 6 USE comvert_mod, ONLY: ap,bp,aps,bps,pseudoalt 7 7 USE logic_mod, ONLY: fxyhypb,ysinus -
trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F
r1422 r1524 3 3 use radinc_h, only: L_NSPECTI 4 4 use radcommon_h, only: WNOI,DWNI 5 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, daysec, dtphys,6 & pi5 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi 6 use time_phylmdz_mod, only: daysec, dtphys 7 7 USE logic_mod, ONLY: fxyhypb,ysinus 8 8 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy -
trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F
r1422 r1524 4 4 use radcommon_h, only: WNOV,DWNV 5 5 use comsoil_h 6 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, daysec, dtphys,7 & pi6 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi 7 use time_phylmdz_mod, only: daysec, dtphys 8 8 USE logic_mod, ONLY: fxyhypb,ysinus 9 9 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy -
trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90
r1482 r1524 19 19 use planete_mod, only: year_day, periastr, apoastr, peri_day, & 20 20 obliquit, z0, lmixmin, emin_turb 21 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, daysec 21 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp 22 use time_phylmdz_mod, only: daysec 22 23 23 24 implicit none -
trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
r1498 r1524 34 34 use planete_mod, only: apoastr, periastr, year_day, peri_day, & 35 35 obliquit, nres, z0 36 use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp, daysec 36 use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp 37 use time_phylmdz_mod, only: daysec 37 38 use callkeys_mod 38 39 implicit none … … 127 128 ! Authors 128 129 ! ------- 129 ! Frederic Hourdin 130 ! Francois Forget 131 ! Christophe Hourdin 130 ! Frederic Hourdin 15/10/93 131 ! Francois Forget 1994 132 ! Christophe Hourdin 02/1997 132 133 ! Subroutine completely rewritten by F. Forget (01/2000) 133 134 ! Water ice clouds: Franck Montmessin (update 06/2003) … … 249 250 250 251 !$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,& 251 252 252 253 !$OMP zdtlw,zdtsw,sensibFlux) 253 254 … … 270 271 real dtlscale(ngrid,nlayer) ! Largescale routine. 271 272 real zdtc(ngrid,nlayer) ! Condense_co2 routine. 272 real zdtdif(ngrid,nlayer) 273 real zdtdif(ngrid,nlayer) ! Turbdiff/vdifc routines. 273 274 real zdtmr(ngrid,nlayer) ! Mass_redistribution routine. 274 275 real zdtrain(ngrid,nlayer) ! Rain routine. … … 442 443 ALLOCATE(tsoil(ngrid,nsoilmx)) 443 444 ALLOCATE(albedo(ngrid,L_NSPECTV)) 444 ALLOCATE(albedo_equivalent(ngrid))445 446 447 445 ALLOCATE(albedo_equivalent(ngrid)) 446 ALLOCATE(albedo_snow_SPECTV(L_NSPECTV)) 447 ALLOCATE(albedo_co2_ice_SPECTV(L_NSPECTV)) 448 ALLOCATE(albedo_bareground(ngrid)) 448 449 ALLOCATE(rnat(ngrid)) 449 450 ALLOCATE(emis(ngrid)) … … 534 535 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 535 536 albedo(:,:)=0.0 536 537 538 537 albedo_bareground(:)=0.0 538 albedo_snow_SPECTV(:)=0.0 539 albedo_co2_ice_SPECTV(:)=0.0 539 540 call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 540 541 … … 736 737 zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi) 737 738 elseif (diurnal) then 738 739 zlss=-2.*pi*(zday-.5) 739 740 else if(diurnal .eqv. .false.) then 740 741 zlss=9999. … … 782 783 zh(ig,l)=pt(ig,l)/zpopsk(ig,l) 783 784 mass(ig,l) = (pplev(ig,l) - pplev(ig,l+1))/glat(ig) 784 785 massarea(ig,l)=mass(ig,l)*area(ig) 785 786 enddo 786 787 enddo … … 845 846 call abort 846 847 endif 847 848 if(water) then 848 849 muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap)) 849 850 851 850 muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) 851 ! take into account water effect on mean molecular weight 852 else 852 853 muvar(1:ngrid,1:nlayer+1)=mugaz 853 854 endif 854 855 855 856 … … 904 905 zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer) 905 906 906 OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV) 907 OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI) 907 OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV) 908 OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI) 908 909 enddo 909 910 … … 978 979 !---------------------------- 979 980 if(enertest)then 980 981 982 981 call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW) 982 call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW) 983 !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk 983 984 call planetwide_sumval(fluxsurfabs_sw(:)*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk 984 985 986 987 985 call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:)/totarea_planet,dEtotsLW) 986 dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:) 987 dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:) 988 if (is_master) then 988 989 print*,'---------------------------------------------------------------' 989 990 print*,'In corrk SW atmospheric heating =',dEtotSW,' W m-2' … … 993 994 print*,'In corrk LW surface heating =',dEtotsLW,' W m-2' 994 995 print*,'surface net rad heating (SW+LW) =',dEtotsLW+dEtotsSW,' W m-2' 995 996 endif 996 997 endif ! end of 'enertest' 997 998 … … 1009 1010 1010 1011 ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff. 1011 1012 1012 if (UseTurbDiff) then 1013 1013 1014 call turbdiff(ngrid,nlayer,nq,rnat, & 1014 1015 ptimestep,capcal,lwrite, & … … 1017 1018 pdt,pdq,zflubid, & 1018 1019 zdudif,zdvdif,zdtdif,zdtsdif, & 1019 1020 sensibFlux,q2,zdqdif,zdqevap,zdqsdif, & 1020 1021 taux,tauy,lastcall) 1021 1022 1022 1023 1023 else 1024 1024 1025 zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer) 1025 1026 … … 1030 1031 zdh,pdq,zflubid, & 1031 1032 zdudif,zdvdif,zdhdif,zdtsdif, & 1032 1033 sensibFlux,q2,zdqdif,zdqsdif, & 1033 1034 taux,tauy,lastcall) 1034 1035 1035 1036 zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only 1036 1037 zdqevap(1:ngrid,1:nlayer)=0. 1037 1038 1038 1039 end if !end of 'UseTurbDiff' … … 1060 1061 if(enertest)then 1061 1062 1062 1063 dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:) 1063 1064 do ig = 1, ngrid 1064 1065 1065 dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground 1066 dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground 1066 1067 enddo 1067 1068 1068 1069 1070 1071 1069 call planetwide_sumval(dEdiff(:)*area(:)/totarea_planet,dEtot) 1070 dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:) 1071 call planetwide_sumval(dEdiffs(:)*area(:)/totarea_planet,dEtots) 1072 call planetwide_sumval(sensibFlux(:)*area(:)/totarea_planet,AtmToSurf_TurbFlux) 1072 1073 1073 1074 if (is_master) then 1074 1075 1075 1076 if (UseTurbDiff) then 1076 1077 1077 print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2' 1078 print*,'In TurbDiff non-cons atm nrj change =',dEtot,' W m-2' 1078 1079 print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2' 1079 1080 1081 1082 1083 1084 1080 else 1081 print*,'In vdifc sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2' 1082 print*,'In vdifc non-cons atm nrj change =',dEtot,' W m-2' 1083 print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2' 1084 end if 1085 endif ! end of 'is_master' 1085 1086 1086 1087 ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere. … … 1091 1092 if(watertest.and.water)then 1092 1093 1093 1094 1094 call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp) 1095 call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots_tmp) 1095 1096 do ig = 1, ngrid 1096 1097 1098 1099 1100 1101 1097 vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap)) 1098 enddo 1099 call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) 1100 call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots) 1101 dWtot = dWtot + dWtot_tmp 1102 dWtots = dWtots + dWtots_tmp 1102 1103 do ig = 1, ngrid 1103 1104 enddo1105 1106 1107 1104 vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice)) 1105 enddo 1106 call planetwide_maxval(vdifcncons(:),nconsMAX) 1107 1108 if (is_master) then 1108 1109 print*,'---------------------------------------------------------------' 1109 1110 print*,'In difv atmospheric water change =',dWtot,' kg m-2' … … 1111 1112 print*,'In difv non-cons factor =',dWtot+dWtots,' kg m-2' 1112 1113 print*,'In difv MAX non-cons factor =',nconsMAX,' kg m-2 s-1' 1113 1114 endif 1114 1115 1115 1116 endif ! end of 'watertest' … … 1157 1158 ! Test energy conservation 1158 1159 if(enertest)then 1159 1160 call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot) 1160 1161 if (is_master) print*,'In convadj atmospheric energy change =',dEtot,' W m-2' 1161 1162 endif … … 1163 1164 ! Test water conservation 1164 1165 if(watertest)then 1165 1166 call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp) 1166 1167 do ig = 1, ngrid 1167 1168 1169 1170 1168 cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap)) 1169 enddo 1170 call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) 1171 dWtot = dWtot + dWtot_tmp 1171 1172 do ig = 1, ngrid 1172 1173 enddo1174 1173 cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice)) 1174 enddo 1175 call planetwide_maxval(cadjncons(:),nconsMAX) 1175 1176 1176 1177 if (is_master) then 1177 1178 print*,'In convadj atmospheric water change =',dWtot,' kg m-2' 1178 1179 print*,'In convadj MAX non-cons factor =',nconsMAX,' kg m-2 s-1' 1179 1180 endif 1180 1181 1181 1182 endif ! end of 'watertest' … … 1208 1209 ! test energy conservation 1209 1210 if(enertest)then 1210 1211 1212 1213 1211 call planetwide_sumval(cpp*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot) 1212 call planetwide_sumval(capcal(:)*zdtsurfc(:)*area(:)/totarea_planet,dEtots) 1213 if (is_master) then 1214 print*,'In co2cloud atmospheric energy change =',dEtot,' W m-2' 1214 1215 print*,'In co2cloud surface energy change =',dEtots,' W m-2' 1215 1216 endif 1216 1217 endif 1217 1218 … … 1233 1234 if(watercond.and.(RLVTT.gt.1.e-8))then 1234 1235 1235 1236 1236 dqmoist(1:ngrid,1:nlayer,1:nq)=0. 1237 dtmoist(1:ngrid,1:nlayer)=0. 1237 1238 1238 1239 ! Moist Convective Adjustment. 1239 1240 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1240 1241 call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man) 1241 1242 1242 1243 pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) & 1243 1244 + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) 1244 1245 pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) & 1245 1246 + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice) … … 1249 1250 if(enertest)then 1250 1251 1251 1252 1253 1254 1252 call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot) 1253 call planetwide_maxval(dtmoist(:,:),dtmoist_max) 1254 call planetwide_minval(dtmoist(:,:),dtmoist_min) 1255 madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:) 1255 1256 do ig=1,ngrid 1256 1257 madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:)) 1257 1258 enddo 1258 1259 1259 1260 if (is_master) then 1260 1261 print*,'In moistadj atmospheric energy change =',dEtot,' W m-2' 1261 1262 print*,'In moistadj MAX atmospheric energy change =',dtmoist_max*ptimestep,'K/step' 1262 1263 print*,'In moistadj MIN atmospheric energy change =',dtmoist_min*ptimestep,'K/step' 1263 1264 1265 call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+&1266 1267 1264 endif 1265 1266 call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+ & 1267 massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) 1268 if (is_master) print*,'In moistadj atmospheric water change =',dWtot,' kg m-2' 1268 1269 1269 1270 endif ! end of 'enertest' 1270 1271 1271 1272 1272 1273 ! Large scale condensation/evaporation. … … 1284 1285 lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:)) 1285 1286 enddo 1286 1287 1288 1287 call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot) 1288 1289 if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2' 1289 1290 1290 1291 ! Test water conservation. 1291 call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+&1292 1292 call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+ & 1293 SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot) 1293 1294 1294 1295 if (is_master) print*,'In largescale atmospheric water change =',dWtot,' kg m-2' 1295 1296 endif ! end of 'enertest' 1296 1297 … … 1316 1317 1317 1318 pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) & 1318 1319 + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) 1319 1320 pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) & 1320 1321 + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice) 1321 1322 pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer) 1322 1323 … … 1327 1328 if(enertest)then 1328 1329 1329 1330 1331 1332 1333 1334 1335 1336 1337 1330 call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot) 1331 if (is_master) print*,'In rain atmospheric T energy change =',dEtot,' W m-2' 1332 call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)/totarea_planet*RLVTT/cpp,dItot_tmp) 1333 call planetwide_sumval(area(:)*zdqssnow(:)/totarea_planet*RLVTT/cpp,dItot) 1334 dItot = dItot + dItot_tmp 1335 call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dVtot_tmp) 1336 call planetwide_sumval(area(:)*zdqsrain(:)/totarea_planet*RLVTT/cpp,dVtot) 1337 dVtot = dVtot + dVtot_tmp 1338 dEtot = dItot + dVtot 1338 1339 1339 1340 if (is_master) then 1340 1341 print*,'In rain dItot =',dItot,' W m-2' 1341 1342 print*,'In rain dVtot =',dVtot,' W m-2' 1342 1343 print*,'In rain atmospheric L energy change =',dEtot,' W m-2' 1343 1344 endif 1344 1345 1345 1346 ! Test water conservation 1346 call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+&1347 1348 1347 call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+ & 1348 massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot) 1349 call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*area(:)*ptimestep/totarea_planet,dWtots) 1349 1350 1350 1351 1352 1353 1354 1351 if (is_master) then 1352 print*,'In rain atmospheric water change =',dWtot,' kg m-2' 1353 print*,'In rain surface water change =',dWtots,' kg m-2' 1354 print*,'In rain non-cons factor =',dWtot+dWtots,' kg m-2' 1355 endif 1355 1356 1356 1357 endif ! end of 'enertest' … … 1373 1374 1374 1375 iq=igcm_h2o_ice 1375 1376 1377 1378 1376 call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot) 1377 call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots) 1378 if (is_master) then 1379 print*,'Before sedim pq =',dWtot,' kg m-2' 1379 1380 print*,'Before sedim pdq =',dWtots,' kg m-2' 1380 1381 endif 1381 1382 endif 1382 1383 … … 1387 1388 if(watertest)then 1388 1389 iq=igcm_h2o_ice 1389 1390 1391 1392 1393 1394 1390 call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot) 1391 call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots) 1392 if (is_master) then 1393 print*,'After sedim pq =',dWtot,' kg m-2' 1394 print*,'After sedim pdq =',dWtots,' kg m-2' 1395 endif 1395 1396 endif 1396 1397 … … 1401 1402 ! Test water conservation 1402 1403 if(watertest)then 1403 1404 1405 1406 1407 1408 1409 1404 call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot) 1405 call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:)*ptimestep/totarea_planet,dWtots) 1406 if (is_master) then 1407 print*,'In sedim atmospheric ice change =',dWtot,' kg m-2' 1408 print*,'In sedim surface ice change =',dWtots,' kg m-2' 1409 print*,'In sedim non-cons factor =',dWtot+dWtots,' kg m-2' 1410 endif 1410 1411 endif 1411 1412 … … 1421 1422 1422 1423 zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) * & 1423 1424 1425 1426 1424 ( zdqevap(1:ngrid,1:nlayer) & 1425 + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap) & 1426 + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap) & 1427 + dqvaplscale(1:ngrid,1:nlayer) ) 1427 1428 1428 1429 do ig = 1, ngrid 1429 1430 zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer)) 1430 1431 enddo 1431 1432 1432 1433 1434 1435 1436 1433 call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr) 1434 call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col) 1435 call writediagfi(ngrid,"mass","mass","kg/m2",3,mass) 1436 1437 call mass_redistribution(ngrid,nlayer,nq,ptimestep, & 1437 1438 rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf, & 1438 1439 1440 1439 pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr, & 1440 zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr) 1441 1441 1442 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq) 1442 1443 dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq) … … 1444 1445 pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer) 1445 1446 pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer) 1446 1447 pdpsrf(1:ngrid) = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid) 1447 1448 zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid) 1448 1449 1449 1450 endif 1450 1451 1451 1452 ! ------------------ … … 1493 1494 call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, & 1494 1495 capcal,albedo,albedo_bareground, & 1495 1496 albedo_snow_SPECTV,albedo_co2_ice_SPECTV, & 1496 1497 mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic, & 1497 1498 sea_ice) … … 1504 1505 ! Test energy conservation 1505 1506 if(enertest)then 1506 1507 1507 call planetwide_sumval(area(:)*capcal(:)*zdtsurf_hyd(:)/totarea_planet,dEtots) 1508 if (is_master) print*,'In hydrol surface energy change =',dEtots,' W m-2' 1508 1509 endif 1509 1510 1510 1511 ! test water conservation 1511 1512 if(watertest)then 1512 1513 1514 1515 1513 call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots) 1514 if (is_master) print*,'In hydrol surface ice change =',dWtots,' kg m-2' 1515 call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots) 1516 if (is_master) then 1516 1517 print*,'In hydrol surface water change =',dWtots,' kg m-2' 1517 1518 print*,'---------------------------------------------------------------' 1518 1519 endif 1519 1520 endif 1520 1521 … … 1590 1591 ! Test energy conservation 1591 1592 if(enertest)then 1592 call planetwide_sumval(area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)1593 1593 call planetwide_sumval(area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots) 1594 if (is_master) print*,'Surface energy change =',dEtots,' W m-2' 1594 1595 endif 1595 1596 … … 1638 1639 print*,Ts1,Ts2,Ts3,TsS 1639 1640 else 1640 1641 if (is_master) then 1641 1642 print*,' ave[Tsurf] min[Tsurf] max[Tsurf]' 1642 1643 print*,Ts1,Ts2,Ts3 1643 1644 endif 1644 1645 end if 1645 1646 … … 1648 1649 if(corrk)then 1649 1650 1650 1651 1652 1653 1654 1651 call planetwide_sumval(area(:)*fluxtop_dn(:)/totarea_planet,ISR) 1652 call planetwide_sumval(area(:)*fluxabs_sw(:)/totarea_planet,ASR) 1653 call planetwide_sumval(area(:)*fluxtop_lw(:)/totarea_planet,OLR) 1654 call planetwide_sumval(area(:)*fluxgrd(:)/totarea_planet,GND) 1655 call planetwide_sumval(area(:)*fluxdyn(:)/totarea_planet,DYN) 1655 1656 do ig=1,ngrid 1656 1657 if(fluxtop_dn(ig).lt.0.0)then … … 1668 1669 DYN=0.0 1669 1670 endif 1670 1671 1671 1672 if (is_master) then 1672 1673 print*,' ISR ASR OLR GND DYN [W m^-2]' 1673 1674 print*, ISR,ASR,OLR,GND,DYN 1674 1675 endif 1675 1676 1676 1677 if(enertest .and. is_master)then … … 1688 1689 open(93,file="tem_bal.out",form='formatted',position='append') 1689 1690 if(callsoil)then 1690 1691 1692 1693 1691 write(93,*) zday,Ts1,Ts2,Ts3,TsS 1692 else 1693 write(93,*) zday,Ts1,Ts2,Ts3 1694 endif 1694 1695 close(93) 1695 1696 endif … … 1744 1745 if(water)then 1745 1746 1746 1747 1748 1749 1747 call planetwide_sumval(area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf) 1748 call planetwide_sumval(area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf) 1749 call planetwide_sumval(area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol) 1750 call planetwide_sumval(area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol) 1750 1751 1751 1752 h2otot = icesrf + liqsrf + icecol + vapcol 1752 1753 1753 1754 if (is_master) then 1754 1755 print*,' Total water amount [kg m^-2]: ',h2otot 1755 1756 print*,' Surface ice Surface liq. Atmos. con. Atmos. vap. [kg m^-2] ' 1756 1757 print*, icesrf,liqsrf,icecol,vapcol 1757 1758 endif 1758 1759 1759 1760 if(meanOLR .and. is_master)then … … 2013 2014 2014 2015 endif ! end of 'callrad' 2015 2016 2016 2017 if(enertest) then 2017 2018 2018 2019 2020 2019 if (calldifv) then 2020 2021 call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2) 2021 2022 call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux) 2022 2023 2023 ! 2024 ! 2025 ! 2026 2027 2024 ! call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff) 2025 ! call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff) 2026 ! call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs) 2027 2028 endif 2028 2029 2029 2030 2031 2032 2030 if (corrk) then 2031 call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw) 2032 call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw) 2033 endif 2033 2034 2034 2035 if(watercond) then 2035 2036 2036 2037 2038 2037 call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 2038 call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE) 2039 call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat) 2039 2040 2040 ! 2041 ! 2042 ! 2041 ! call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) 2042 ! call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz) 2043 ! call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol) 2043 2044 2044 2045 endif … … 2060 2061 ! Output aerosols. 2061 2062 if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) & 2062 2063 call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2)) 2063 2064 if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & 2064 2065 call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o)) 2065 2066 if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) & 2066 2067 call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2)) 2067 2068 if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) & 2068 2069 call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o)) 2069 2070 2070 2071 ! Output tracers. … … 2086 2087 call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac) 2087 2088 call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac) 2088 2089 call writediagfi(ngrid,"RH","relative humidity"," ",3,RH) 2089 2090 endif 2090 2091 -
trunk/LMDZ.GENERIC/libf/phystd/soil.F
r1516 r1524 5 5 6 6 use comsoil_h, only: layer, mlayer, volcapa, inertiedat 7 use comcstfi_mod, only: pi, daysec 7 use comcstfi_mod, only: pi 8 use time_phylmdz_mod, only: daysec 8 9 use planete_mod, only: year_day 9 10 use comgeomfi_h, only: long, lati -
trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
r1482 r1524 53 53 use planete_mod, only: year_day, periastr, apoastr, peri_day, 54 54 & obliquit, z0, lmixmin, emin_turb 55 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, dtphys,56 & daysec, r55 use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r 56 use time_phylmdz_mod, only: dtphys, daysec 57 57 use callkeys_mod, only: check_cpp_match,force_cpp 58 58 implicit none … … 507 507 write(*,*) 508 508 509 ENDIF !of if (Lmodif == 1)509 ENDIF ! of if (Lmodif == 1) 510 510 511 511 c----------------------------------------------------------------------- -
trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F
r1422 r1524 44 44 & is_master, gather 45 45 USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo 46 USE t emps_mod, ONLY: day_ini46 USE time_phylmdz_mod, ONLY: day_ini 47 47 implicit none 48 48 -
trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F
r1422 r1524 50 50 use control_mod, only: ecritphy, iphysiq, day_step 51 51 use callkeys_mod, only: iradia 52 USE temps_mod, ONLY: day_ini52 use time_phylmdz_mod, ONLY: day_ini 53 53 54 54 implicit none -
trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F
r1422 r1524 50 50 use control_mod, only: ecritphy, iphysiq, day_step 51 51 use callkeys_mod, only: iradia 52 USE temps_mod, ONLY: day_ini52 use time_phylmdz_mod, ONLY: day_ini 53 53 54 54 implicit none -
trunk/LMDZ.GENERIC/libf/phystd/writeg1d.F
r1397 r1524 129 129 130 130 SUBROUTINE endg1d(ngrid,nlayer,zlayer,ndt) 131 USE comcstfi_mod, ONLY: dtphys, daysec131 USE time_phylmdz_mod, ONLY: dtphys, daysec 132 132 USE comg1d_mod, ONLY: g1d_nomfich,g1d_unitfich,g1d_nvar, 133 133 & g1d_nomvar,g1d_titrevar,g1d_dimvar,g1d_nlayer,g1d_unitctl,
Note: See TracChangeset
for help on using the changeset viewer.
