source: trunk/LMDZ.PLUTO.old/libf/phypluto/testphys1d_triton.F @ 3175

Last change on this file since 3175 was 3175, checked in by emillour, 10 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM
File size: 26.1 KB
RevLine 
[3175]1      PROGRAM testphys1d_triton
2
3! to use  'getin'
4      USE ioipsl_getincom
5      use planet_h         
6      use comgeomfi_h
7      use comsoil_h
8      IMPLICIT NONE
9
10!==================================================================
11!     
12!     Purpose
13!     -------
14!     Run the physics package of the universal model in a 1D column.
15!     
16!     It can be compiled with a command like (e.g. for 25 layers):
17!                                  "makegcm -p pluto -d 25 testphys1d"
18!     It requires the files "callphys.def", "z2sig.def",
19!     "traceur.def", and "run1d.def" with a line "INCLUDEDEF=callphys.def"
20!
21!     Authors
22!     -------
23!     Frederic Hourdin
24!     R. Fournier
25!     F. Forget
26!     F. Montmessin (water ice added)
27
28!==================================================================
29
30#include "dimensions.h"
31#include "dimphys.h"
32#include "surfdat.h"
33!#include "comsoil.h"
34#include "comdiurn.h"
35#include "callkeys.h"
36#include "comcstfi.h"
37#include "comsaison.h"
38#include "control.h"
39#include "comvert.h"
40#include "netcdf.inc"
41#include "comg1d.h"
42#include "fisice.h"
43#include "logic.h"
44#include "advtrac.h"
45
46c --------------------------------------------------------------
47c  Declarations
48c --------------------------------------------------------------
49c
50      INTEGER nlayer,nlevel,nsoil,ndt
51      INTEGER ilayer,ilevel,isoil,idt,iq
52      LOGICAl firstcall,lastcall
53c
54      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
55      INTEGER lecttsoil     ! lecture of tsoil from proftsoil
56      INTEGER lecthaze      ! lecture of haze from profhaze
57      REAL day           ! date durant le run
58      REAL AU           ! astronomical unit AU=149 million km
59      REAL time             ! time (0<time<1 ; time=0.5 a midi)
60      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
61      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
62      REAL psurf,tsurf     
63      REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
64      REAL gru,grv   ! prescribed "geostrophic" background wind
65      REAL temp(nlayermx)   ! temperature at the middle of the layers
66      REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
67      REAL qsurf(nqmx)      ! tracer surface budget (e.g. kg.m-2)
68      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
69      integer :: i_n2=0  ! tracer index of n2 ice
70      integer :: i_ch4_ice=0  ! tracer index of Ch4 ice
71      integer :: i_ch4_gas=0  ! tracer index of Ch4 gas
72      integer :: i_co_ice=0  ! tracer index of Co ice
73      integer :: i_co_gas=0  ! tracer index of Co gas
74      integer :: i_prec_haze=0  ! tracer index of haze precursor
75      integer :: i_haze=0  ! tracer index of haze
76      integer :: i_haze10=0  ! tracer index of haze
77      integer :: i_haze30=0  ! tracer index of haze
78      integer :: i_haze50=0  ! tracer index of haze
79      integer :: i_haze100=0  ! tracer index of haze
80      REAL emis             ! surface layer
81      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
82      REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
83      REAL  Lsperi,excentric
84
85c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
86      REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
87      REAL dpsurf   
88      REAL dq(nlayermx,nqmx)
89
90c   Various intermediate variables
91      REAL zls
92      REAL phi(nlayermx),h(nlayermx),s(nlayermx)
93      REAL pks, ptif, w(nlayermx)
94      INTEGER ierr
95      REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
96      Logical  tracerdyn
97      integer :: nq=1 ! number of tracers
98
99      character (len=7) :: str7
100
101      logical saveprofile
102
103!     added by RW for zlay computation
104      real Hscale, rho, dz
105
106      real ch4ref ! % ch4
107      real coref ! % co
108      real hazeref ! haze kg/kg
109      real prechazeref ! prec haze kg/kg
110
111
112c=======================================================================
113
114c=======================================================================
115c INITIALISATION
116c=======================================================================
117
118      saveprofile=.true.
119
120c ------------------------------------------------------
121c  Default values for constants (corresponding to Triton)
122c ------------------------------------------------------
123
124      pi=2.E+0*asin(1.E+0)
125
126c     Constante de Triton
127c     ----------------------------
128      rad=1353400.               ! rayon de Triton (m)
129      daysec=507773            ! duree du sol (s)
130      omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
131      g=0.78                 ! gravite (m.s-2)
132      mugaz=28.                ! Masse molaire de l'atm (g.mol-1) ~28
133      rcp=.2853                ! = r/cp  ~0.2853
134      r= 8.314511E+0 *1000.E+0/mugaz
135      cpp= r/rcp
136      year_day = 10247         ! duree de l'annee (sols)
137      periheli = 4444.45       ! dist.min. soleil-Triton (Mkm)
138      aphelie = 4545.67         ! dist.max. soleil-Triton (Mkm)
139      peri_day =  0      ! date du perihelie (sols depuis printemps)
140      obliquit=28.32         ! Obliquite de la planete (deg)         
141      excentric= 0.      ! Excentricite
142 
143c     Parametres Couche limite et Turbulence
144c     --------------------------------------
145      z0 =  1.e-2                ! surface roughness (m) ~0.01
146 
147c     ----------------------------------------------------
148      hybrid=.true.
149
150!     Constants for stokes.F90
151      molrad=1.93e-10 ! N2 TB16!!!
152      visc=6.67e-6     ! N2 TB16!!!
153
154c ------------------------------------------------------
155c  Load parameters from "run.def"
156c ------------------------------------------------------
157      open(90,file='run.def',status='old',form='formatted',
158     &     iostat=ierr)
159      if (ierr.ne.0) then
160        write(*,*) 'Cannot find required file "run.def"'
161        write(*,*) '  (which should contain some input parameters'
162        write(*,*) '   along with the following line:'
163        write(*,*) '   INCLUDEDEF=callphys.def'
164        write(*,*) '   )'
165        write(*,*) ' ... might as well stop here ...'
166        stop
167      else
168        close(90)
169      endif
170
171! check if we are going to run with or without tracers
172      write(*,*) "Run with or without tracer transport ?"
173      tracer=.false. ! default value
174      call getin("tracer",tracer)
175      write(*,*) " tracer = ",tracer
176
177! while we're at it, check if there is a 'traceur.def' file
178! and preocess it, if necessary. Otherwise initialize tracer names
179      if (tracer) then
180      ! load tracer names from file 'traceur.def'
181        open(90,file='traceur.def',status='old',form='formatted',
182     &       iostat=ierr)
183        if (ierr.eq.0) then
184          write(*,*) "testphys1d: Reading file traceur.def"
185          ! read number of tracers:
186          read(90,*,iostat=ierr) nq
187          if (ierr.ne.0) then
188            write(*,*) "testphys1d: error reading number of tracers"
189            write(*,*) "   (first line of traceur.def) "
190            stop
191          else
192            ! check that the number of tracers is indeed nqmx
193            if (nq.ne.nqmx) then
194              write(*,*) "testphys1d: error, wrong number of tracers:"
195              write(*,*) "nq=",nq," whereas nqmx=",nqmx
196              stop
197            endif
198          endif
199       
200          ! initialize advection schemes to Van-Leer for all tracers
201          do iq=1,nq
202            iadv(iq)=3 ! Van-Leer
203          enddo
204       
205          do iq=1,nq
206          ! minimal version, just read in the tracer names, 1 per line
207            read(90,*,iostat=ierr) tnom(iq)
208            write(*,*) 'Ini 1d reading traceur.def: ', tnom(iq),iq
209            if (ierr.ne.0) then
210              write(*,*) 'testphys1d: error reading tracer names...'
211              stop
212            endif
213          enddo !of do iq=1,nq
214
215! TB check for n2_ice / ch4 tracers:
216         write(*,*) 'Tracers: tnom=', tnom(:)
217         do iq=1,nq
218           if (tnom(iq)=="n2") then
219             i_n2=iq
220           elseif (tnom(iq)=="ch4_ice") then
221             i_ch4_ice=iq
222           elseif (tnom(iq)=="ch4_gas") then
223             i_ch4_gas=iq
224           elseif (tnom(iq)=="prec_haze") then
225             i_prec_haze=iq
226           elseif (tnom(iq)=="haze") then
227             i_haze=iq
228           elseif (tnom(iq)=="haze10") then
229             i_haze10=iq
230           elseif (tnom(iq)=="haze30") then
231             i_haze30=iq
232           elseif (tnom(iq)=="haze50") then
233             i_haze50=iq
234           elseif (tnom(iq)=="haze100") then
235             i_haze100=iq
236           elseif (tnom(iq)=="co_ice") then
237             i_co_ice=iq
238           elseif (tnom(iq)=="co_gas") then
239             i_co_gas=iq
240           endif
241         enddo
242
243        else  ! ierr=0
244          write(*,*) 'Cannot find required file "traceur.def"'
245          write(*,*) ' If you want to run with tracers, I need it'
246          write(*,*) ' ... might as well stop here ...'
247          stop
248        endif
249        close(90)
250
251      else
252      ! we still need to set (dummy) tracer names for physdem1
253        nq=nqmx
254        do iq=1,nq
255          write(str7,'(a1,i2.2)')'q',iq
256          tnom(iq)=str7
257        enddo
258        ! actually, we'll need at least one "n2_ice" tracer
259        ! (for surface N2 ice)
260        write(*,*) "No tracer ! tracer=false"
261        tnom(1)="n2"
262        i_n2=1
263      endif ! of if (tracer)
264
265c  Date et heure locale du debut du run
266c  ------------------------------------
267c  Date (en sols depuis le solstice de printemps) du debut du run
268      day0 = 0 ! default value for day0
269      write(*,*) 'Initial date (in sols ; =0 at Ls=0)?'
270      call getin("day0",day0)
271      day=float(day0)
272      write(*,*) " day0 = ",day0
273c  Heure de demarrage
274      time=0 ! default value for time
275      write(*,*)'Initial local time (in hours, between 0 and 24)?'
276      call getin("time",time)
277      write(*,*)" time = ",time
278      time=time/24.E+0 ! convert time (hours) to fraction of sol
279
280c  Discretisation (Definition de la grille et des pas de temps)
281c  --------------
282c
283      nlayer=nlayermx
284      nlevel=nlayer+1
285      nsoil=nsoilmx
286      PRINT *,'Dims nlevel=',nlevel,' nsoil=',nsoil
287
288      PRINT *,'Length of day (s) ?'
289      call getin("daysec",daysec)
290      write(*,*) " daysec = ",daysec
291      omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
292
293      PRINT *,'Length of year (sol) ?'
294      call getin("year_day",year_day)
295      write(*,*) " year_day= ",year_day
296
297      day_step=48 ! default value for day_step
298      PRINT *,'Number of time steps per sol ?'
299      call getin("day_step",day_step)
300      write(*,*) " day_step = ",day_step
301
302      ndt=10 ! default value for ndt
303      PRINT *,'Number of sols to run ?'
304      call getin("ndt",ndt)
305      write(*,*) " ndt = ",ndt
306
307      ndt=ndt*day_step     
308      dtphys=daysec/day_step 
309
310c Pression de surface sur la planete
311c ------------------------------------
312      psurf=1. ! default value for psurf
313      PRINT *,'Surface pressure (Pa) ?'
314      call getin("psurf",psurf)
315      write(*,*) " psurf = ",psurf
316c Pression de reference
317      preff=1. ! these values are not needed in 1D 
318      pa=0.25*preff
319c Gravity of planet
320c -----------------
321      g=0.78 ! default value for g
322      PRINT *,'Gravity ?'
323      call getin("g",g)
324      write(*,*) " g = ",g
325
326c Molar mass of atmosphere
327c ------------------------
328      PRINT *,'Molar mass of atmosphere ?'
329      call getin("mugaz",mugaz)
330      write(*,*) " mugaz = ",mugaz
331
332c Specific heat capacity of atmosphere
333c ------------------------------------
334      PRINT *,'Specific heat capacity of atmosphere ?'
335      call getin("cpp",cpp)
336      write(*,*) " cpp = ",cpp
337
338      r= 8.314511E+0 *1000.E+0/mugaz
339      rcp=r/cpp
340 
341c  latitude/longitude
342c  ------------------
343      lati(1)=0 ! default value for lati(1)
344      PRINT *,'latitude (in degrees) ?'
345      call getin("latitude",lati(1))
346      write(*,*) " latitude = ",lati(1)
347      lati(1)=lati(1)*pi/180.E+0
348      long(1)=0.E+0
349      long(1)=long(1)*pi/180.E+0
350
351c  periastron/apoastron
352c  --------------------
353      PRINT *,'periastron (in AU) ?'
354      call getin("periheli",periheli)
355      write(*,*) "periastron = ",periheli
356      AU=149.597870700 ! km
357      periheli=periheli*AU ! AU to Mkm
358      PRINT *,'apoastron (in AU) ?'
359      call getin("aphelie",aphelie)
360      write(*,*) "apoastron = ",aphelie
361      aphelie=aphelie*AU ! AU to Mkm
362
363      excentric=(1-periheli*periheli/(aphelie*aphelie) )**0.5
364
365c  obliquity
366c  ---------
367      PRINT *,'obliquity (in deg) ?'
368      call getin("obliquit",obliquit)
369      write(*,*) "obliquity = ",obliquit
370
371c  Initialisation albedo /emis / inertie du sol
372c  --------------------------------------
373      albedodat(1)=0.2 ! default value for albedodat
374      PRINT *,'Albedo of bare ground ?'
375      call getin("albedo",albedodat(1))
376      write(*,*) " albedo = ",albedodat(1)
377
378      emis=0.8 ! default value for emissivity
379      PRINT *,'Emissivity of bare ground ?'
380      call getin("emis",emis)
381      write(*,*) " emis = ",emis
382
383      inertiedat(1,1)=400 ! default value for inertiedat
384      PRINT *,'Soil thermal inertia (SI) ?'
385      call getin("inertia",inertiedat(1,1))
386      write(*,*) " inertia = ",inertiedat(1,1)
387c
388c  pour le schema d'ondes de gravite
389c  ---------------------------------
390      zmea(1)=0.E+0
391      zstd(1)=0.E+0
392      zsig(1)=0.E+0
393      zgam(1)=0.E+0
394      zthe(1)=0.E+0
395
396c    Initialisation des traceurs
397c    ---------------------------
398      DO iq=1,nqmx
399        if (iq.eq.i_n2) then
400           DO ilayer=1,nlayer
401              q(ilayer,iq) = 1
402           ENDDO
403           write(*,*) 'ini 1D traceur ',iq,' : q_n2 = ', q(:,iq)
404        else if (iq.eq.i_ch4_gas) then
405          ch4ref=0.5 ! default value for vmr ch4
406          PRINT *,'vmr CH4 (%) ?'
407          call getin("ch4ref",ch4ref)
408          write(*,*) " CH4 ref (%) = ",ch4ref
409          DO ilayer=1,nlayer
410             q(ilayer,iq) = 0.01*ch4ref*16./28.
411          ENDDO
412          write(*,*) 'ini 1D traceur ',iq,' : q_ch4_gas = ', q(:,iq)
413        else if (iq.eq.i_co_gas) then
414          coref=0.05 ! default value for vmr ch4
415          PRINT *,'vmr CO (%) ?'
416          call getin("coref",coref)
417          write(*,*) " CO ref (%) = ",coref
418          DO ilayer=1,nlayer
419             q(ilayer,iq) = 0.01*coref*28./28.
420          ENDDO
421          write(*,*) 'ini 1D traceur ',iq,' : q_co_gas = ', q(:,iq)
422        else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30).or.
423     &     (iq.eq.i_haze).or.(iq.eq.i_haze50).or.(iq.eq.i_haze100)) then
424          hazeref=0. ! default value for haze mmr
425          PRINT *,'qhaze (kg/kg) ?'
426          call getin("hazeref",hazeref)
427          write(*,*) " haze ref (kg/kg) = ",hazeref
428          DO ilayer=1,nlayer
429             q(ilayer,iq) = hazeref
430          ENDDO
431          write(*,*) 'ini 1D traceur ',iq,' : q_haze = ', q(:,iq)
432        else if (iq.eq.i_prec_haze) then
433          prechazeref=0. ! default value for vmr ch4
434          PRINT *,'qprechaze (kg/kg) ?'
435          call getin("prechazeref",prechazeref)
436          write(*,*) " prec haze ref (kg/kg) = ",prechazeref
437          DO ilayer=1,nlayer
438             q(ilayer,iq) = prechazeref
439          ENDDO
440          write(*,*) 'ini 1D traceur ',iq,' : q_prec_haze = ', q(:,iq)
441         
442        else
443           DO ilayer=1,nlayer
444              q(ilayer,iq) = 0.
445           ENDDO
446        endif
447      ENDDO
448
449c     TB: clean surface
450      DO iq=1,nqmx
451        qsurf(iq) = 0.
452      ENDDO
453
454c   Initialisation speciales "physiq"
455c   ---------------------------------
456c   la surface de chaque maille est inutile en 1D --->
457      area(1)=1.E+0
458
459c   le geopotentiel au sol est inutile en 1D car tout est controle
460c   par la pression de surface --->
461      phisfi(1)=0.E+0
462
463c   "inifis" reproduit un certain nombre d'initialisations deja faites
464c  + lecture des clefs de callphys.def
465c  + calcul de la frequence d'appel au rayonnement
466c  + calcul des sinus et cosinus des longitude latitude
467
468!   Possible issue with dtphys in input and include!!!
469      CALL inifis(1,llm,day0,daysec,dtphys,
470     .            lati,long,area,rad,g,r,cpp)
471c   Initialisation pour prendre en compte les vents en 1-D
472c   ------------------------------------------------------
473      ptif=2.E+0*omeg*sinlat(1)
474
475c    vent geostrophique
476      gru=10. ! default value for gru
477      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
478      call getin("u",gru)
479      write(*,*) " u = ",gru
480      grv=0. !default value for grv
481      PRINT *,'meridional northward component of the geostrophic',
482     &' wind (m/s) ?'
483      call getin("v",grv)
484      write(*,*) " v = ",grv
485
486c     Initialisation des vents  au premier pas de temps
487      DO ilayer=1,nlayer
488         u(ilayer)=gru
489         v(ilayer)=grv
490      ENDDO
491
492c     energie cinetique turbulente
493      DO ilevel=1,nlevel
494         q2(ilevel)=0.E+0
495      ENDDO
496
497c     Surface
498c     -------------------------------------------
499      if(i_n2.gt.0)then
500         qsurf(i_n2)=0 ! default value for N2ice
501         print*,'Initial N2 ice on the surface (kg.m-2)'
502         call getin("n2ice",qsurf(i_n2))
503         write(*,*) " n2ice = ",qsurf(i_n2)
504      endif
505
506c  calcul des pressions et altitudes en utilisant les niveaux sigma
507c  ----------------------------------------------------------------
508
509c    Vertical Coordinates
510c    """"""""""""""""""""
511      hybrid=.false.
512      PRINT *,'Hybrid coordinates ?'
513      call getin("hybrid",hybrid)
514      write(*,*) " hybrid = ", hybrid
515
516      CALL disvert
517
518      ! we want only the scale height from z2sig, in order to compute zlay
519      open(91,file="z2sig.def",form='formatted')
520      read(91,*) Hscale
521      close(91)
522
523      DO ilevel=1,nlevel
524        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
525      ENDDO
526
527      DO ilayer=1,nlayer
528        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
529      ENDDO
530
531      DO ilayer=1,nlayer
532!        zlay(ilayer)=-300.E+0 *r*log(play(ilayer)/plev(1))
533!     &   /g
534        zlay(ilayer)=-1000.0*Hscale*log(play(ilayer)/plev(1))
535      ENDDO
536
537c   Orbital parameters
538c   ------------------
539      PRINT *,'Solar longitude of perihelion Lsperi '
540      call getin("Lsperi",Lsperi)
541      write(*,*) " Lsperi = ", Lsperi
542     
543      Lsperi = Lsperi * pi/180.   ! Ls of perihelion
544
545c     Calcul du sol correspondant a Lsperi
546      call call_dayperi(Lsperi,excentric,peri_day,year_day)
547      PRINT *,'date of perihelion (sol)',peri_day
548
549c   Profil de temperature au premier appel
550c   --------------------------------------
551      pks=psurf**rcp
552
553c   Altitude en km dans profile: on divise zlay par 1000
554      tmp1(0)=0.E+0
555      DO ilayer=1,nlayer
556        tmp1(ilayer)=zlay(ilayer)/1000.E+0
557      ENDDO
558      call profile(nlayer+1,tmp1,tmp2)
559
560      tsurf=tmp2(0)
561      DO ilayer=1,nlayer
562        temp(ilayer)=tmp2(ilayer)
563      ENDDO
564
565!   Initialize soil properties and temperature
566!   ------------------------------------------
567      volcapa=1.e6 ! volumetric heat capacity
568
569      lecttsoil=0 ! default value for lecttsoil
570      call getin("lecttsoil",lecttsoil)
571      if (lecttsoil==1) then
572         OPEN(14,file='proftsoil',status='old',form='formatted',err=101)
573         DO isoil=1,nsoil
574            READ (14,*) tsoil(isoil)
575            inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
576         ENDDO
577         GOTO 201
578101      STOP'fichier proftsoil inexistant'
579201      CONTINUE
580         CLOSE(14)
581
582      else
583        DO isoil=1,nsoil
584         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
585         tsoil(isoil)=tsurf  ! soil temperature
586        ENDDO
587      endif
588
589!   Initialize depths
590!   -----------------
591      do isoil=0,nsoil-1
592        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
593      enddo
594      do isoil=1,nsoil
595        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
596      enddo
597
598c   Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
599c   ----------------------------------------------------------------
600c   (effectuee avec "writeg1d" a partir de "physiq.F"
601c   ou tout autre subroutine physique, y compris celle ci).
602
603      g1d_nlayer=nlayer
604      g1d_nomfich='g1d.dat'
605      g1d_unitfich=40
606      g1d_nomctl='g1d.ctl'
607      g1d_unitctl=41
608      g1d_premier=.true.
609      g2d_premier=.true.
610
611!     Initialize haze profile
612!     ------------------------------------------
613      if (haze) then
614
615        lecthaze=0 ! default value for lecthaze
616        call getin("lecthaze",lecthaze)
617        if (lecthaze==1) then
618         OPEN(15,file='profhaze',status='old',form='formatted',err=301)
619         DO iq=1,nq
620          if (iq.eq.i_haze) then
621            DO ilayer=1,nlayer
622              READ (15,*) q(ilayer,iq)
623            ENDDO
624          endif
625         ENDDO
626         GOTO 401
627301      STOP'fichier profhaze inexistant'
628401      CONTINUE
629         CLOSE(15)
630        endif
631      endif
632
633c   Ecriture de "startfi"
634c   --------------------
635c   (Ce fichier sera aussitot relu au premier
636c   appel de "physiq", mais il est necessaire pour passer
637c   les variables purement physiques a "physiq"...
638
639      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
640     .              dtphys,float(day0),
641     .              time,tsurf,tsoil,emis,q2,qsurf,
642     .              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
643
644c=======================================================================
645c=======================================================================
646c=======================================================================
647c  BOUCLE TEMPORELLE DU MODELE 1D
648c=======================================================================
649c=======================================================================
650c=======================================================================
651
652      firstcall=.true.
653      lastcall=.false.
654
655      DO idt=1,ndt
656        IF (idt.eq.ndt-day_step-1) then       !test
657         lastcall=.true.
658         call solarlong(day*1.0,zls)
659         write(103,*) 'Ls=',zls*180./pi
660         write(103,*) 'Lat=', lati(1)*180./pi
661         write(103,*) 'RunEnd - Atmos. Temp. File'
662         write(103,*) 'RunEnd - Atmos. Temp. File'
663         write(104,*) 'Ls=',zls*180./pi
664         write(104,*) 'Lat=', lati(1)
665         write(104,*) 'RunEnd - Atmos. Temp. File'
666        ENDIF
667
668c    calcul du geopotentiel
669c     ~~~~~~~~~~~~~~~~~~~~~
670        DO ilayer=1,nlayer
671            s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
672            h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
673        ENDDO
674
675        phi(1)=pks*h(1)*(1.E+0-s(1))
676        DO ilayer=2,nlayer
677          phi(ilayer)=phi(ilayer-1)+
678     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
679     &                  *(s(ilayer-1)-s(ilayer))
680        ENDDO
681
682c       appel de la physique
683c       --------------------
684        CALL physiq (1,llm,nqmx,
685     ,     firstcall,lastcall,
686     ,     day,time,dtphys,
687     ,     plev,play,phi,
688     ,     u, v,temp, q, 
689     ,     w,
690C - sorties
691     s     du, dv, dtemp, dq,dpsurf,tracerdyn)
692
693
694c       evolution du vent : modele 1D
695c       -----------------------------
696 
697c       la physique calcule les derivees temporelles de u et v.
698c       on y rajoute betement un effet Coriolis.
699c
700c       DO ilayer=1,nlayer
701c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
702c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
703c       ENDDO
704
705c       Pour certain test : pas de coriolis a l'equateur
706c       if(lati(1).eq.0.) then
707        DO ilayer=1,nlayer
708             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
709             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
710        ENDDO
711c       end if
712c     
713c       Calcul du temps au pas de temps suivant
714c       ---------------------------------------
715        firstcall=.false.
716        time=time+dtphys/daysec
717        IF (time.gt.1.E+0) then
718            time=time-1.E+0
719            day=day+1
720        ENDIF
721
722c       calcul des vitesses et temperature au pas de temps suivant
723c       ----------------------------------------------------------
724
725        DO ilayer=1,nlayer
726           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
727           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
728           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
729        ENDDO
730
731c       calcul des pressions au pas de temps suivant
732c       ----------------------------------------------------------
733
734        psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
735        DO ilevel=1,nlevel
736          plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
737        ENDDO
738        DO ilayer=1,nlayer
739          play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
740        ENDDO
741
742c       calcul traceur au pas de temps suivant
743c       --------------------------------------
744        DO iq = 1, nqmx
745          DO ilayer=1,nlayer
746             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
747          ENDDO
748        END DO
749
750      ENDDO   ! fin de la boucle temporelle
751
752c    ========================================================
753c    GESTION DES SORTIE
754c    ========================================================
755
756c    fermeture pour conclure les sorties format grads dans "g1d.dat"
757c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
758c    ou tout autre subroutine physique, y compris celle ci).
759
760! save atm temperature profile
761      if(saveprofile)then
762         OPEN(12,file='profile.out',form='formatted')
763         write(12,*) tsurf
764         DO ilayer=1,nlayermx
765            write(12,*) temp(ilayer)
766         ENDDO
767         CLOSE(12)
768      endif
769! save haze profile
770      if (haze.and.lecthaze.eq.1) then
771         OPEN(16,file='profhaze.out',form='formatted')
772         DO iq=1,nq
773          if (iq.eq.i_haze) then
774            DO ilayer=1,nlayer
775              write(16,*) q(ilayer,iq)
776            ENDDO
777          endif
778         ENDDO
779         CLOSE(16)
780      endif
781
782      ! here we compute altitude CORRECTLY!!!
783
784!      print*,'zlay=',zlay/1000.
785!      rho = play(1)/(R*tsurf)
786!      zlay(1)=(plev(1)-plev(2))/(g*rho)
787!      do ilayer=2,nlayer
788!         rho  = play(ilayer)/(R*temp(ilayer))
789!         rho  = play(ilayer)/(R*230.0)
790!         dz = (plev(ilayer-1)-plev(ilayer))/(g*rho)
791!         zlay(ilayer)=zlay(ilayer-1)+dz
792!      enddo
793!      print*,'zlay=',zlay/1000.
794!      stop
795
796c        CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
797      do ilayer=1,nlayer
798        zlay(ilayer)=-18.0e3*log(play(ilayer)/psurf)
799      enddo
800      CALL endg1d(1,nlayer,zlay/1000.,ndt)
801
802c    ========================================================
803      END
804 
805c***********************************************************************
806c***********************************************************************
807c     Subroutines Bidons utilise seulement en 3D, mais
808c     necessaire a la compilation de testphys1d en 1-D
809
810      subroutine gr_fi_dyn
811      RETURN
812      END
813 
814c***********************************************************************
815c***********************************************************************
816
817#include "../dyn3d/disvert.F"
818#include "call_dayperi.F"
Note: See TracBrowser for help on using the repository browser.