source: trunk/LMDZ.VENUS/Tools/lat_torques.F90 @ 816

Last change on this file since 816 was 816, checked in by slebonnois, 12 years ago

SL: tools for postprocessing (Veznus and Titan); see DOC/documentation/vt-tools.pdf
File size: 23.7 KB
Line 
1
2
3program torques
4
5! SL 12/2009:
6! This program reads 4D (lon-lat-alt-time) fields directly from CAM outputs
7! without regrid
8!
9! it computes:
10! dmass -- 4D -- mass of each cell
11! osam  -- 4D -- specific angular momentum (omega term)
12! rsam  -- 4D -- specific angular momentum (zonal wind term)
13! oaam  -- 1D -- integrated angular momentum (omega term)
14! raam  -- 1D -- integrated angular momentum (zonal wind term)
15! tmou  -- 2D -- Mountain torque fonction of lat
16! tbls  -- 2D -- Surface friction torque IF duvdf is present
17!                                       or if simple friction
18! tdyn  -- 2D -- Dynamics torque IF dudyn is present
19!
20! Minimal requirements and dependencies:
21! The dataset must include the following data:
22! - surface pressure and surface geopotential
23! - zonal wind
24! Optional: dudyn, duvdf
25
26implicit none
27
28include "netcdf.inc" ! NetCDF definitions
29
30character (len=128) :: infile ! input file name (name_P.nc)
31character (len=128) :: outfile ! output file name
32
33character (len=64) :: text ! to store some text
34integer infid ! NetCDF input file ID
35integer outfid ! NetCDF output file ID
36integer lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
37integer lon_varid,lat_varid,alt_varid,time_varid, tmpvarid
38integer              :: datashape1d ! shape of 1D datasets
39integer,dimension(2) :: datashape2d ! shape of 2D datasets
40integer,dimension(3) :: datashape3d ! shape of 3D datasets
41integer,dimension(4) :: datashape4d ! shape of 4D datasets
42
43real :: miss_val=9.99e+33 ! special "missing value" to specify missing data
44real,parameter :: miss_val_def=9.99e+33 ! default value for "missing value"
45real :: pi
46real,dimension(:),allocatable :: lon ! longitude
47integer lonlength ! # of grid points along longitude
48real,dimension(:),allocatable :: lat ! latitude
49real,dimension(:),allocatable :: latrad ! latitude in radian
50integer latlength ! # of grid points along latitude
51real,dimension(:),allocatable :: plev ! Pressure levels (Pa)
52integer altlength ! # of grid point along presnivs (of input datasets)
53real,dimension(:),allocatable :: time ! time
54integer timelength ! # of points along time
55real :: p0 ! GCM reference pressure
56real,dimension(:,:,:),allocatable :: ps ! surface pressure
57real,dimension(:,:,:),allocatable :: phis3d ! surface geopotential
58real,dimension(:,:),allocatable :: phis ! surface geopotential
59real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
60real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure
61real,dimension(:),allocatable :: pp ! Pressure levels (Pa)
62real,dimension(:,:,:,:),allocatable :: temp ! atmospheric temperature
63real,dimension(:,:,:,:),allocatable :: vitu ! zonal wind (in m/s)
64real,dimension(:,:,:,:),allocatable :: vitv ! meridional wind (in m/s)
65real,dimension(:,:,:,:),allocatable :: vitw ! vertical wind (in Pa/s, then converted in m/s)
66real,dimension(:,:,:,:),allocatable :: duvdf  ! Friction in BL
67real,dimension(:,:,:,:),allocatable :: dudyn  ! Dynamics
68
69real,dimension(:,:,:,:),allocatable :: grav ! gravity field (m s-2)
70real,dimension(:,:,:,:),allocatable :: dmass ! mass in cell (kg)
71real,dimension(:,:,:,:),allocatable :: osam ! planetary rotation specific ang (m2/s)
72real,dimension(:,:,:,:),allocatable :: rsam ! zonal wind specific ang (m2/s)
73real,dimension(:),allocatable :: oaam ! planetary rotation total ang (kg m2/s)
74real,dimension(:),allocatable :: raam ! zonal wind total ang (kg m2/s)
75real,dimension(:,:),allocatable :: tmou ! mountain torque (kg m2/s2)
76real,dimension(:,:),allocatable :: tdyn ! dynamics torque (kg m2/s2)
77real,dimension(:,:),allocatable :: tbls ! friction torque (kg m2/s2)
78
79! for angular momentum budget normalisation
80real,parameter :: hadley=1.e18
81real,parameter :: hadday=1.e25
82
83integer ierr,ierr1,ierr2 ! NetCDF routines return codes
84integer i,j,ilon,ilat,ilev,itim ! for loops
85integer idlsurf ! for option ideal surface
86logical :: flag_duvdf,flag_dudyn,lmdflag
87
88real :: deltalat,deltalon ! lat and lon intervals in radians
89real,dimension(:,:,:,:),allocatable :: deltap ! pressure thickness of each layer (Pa)
90real :: tmpp ! temporary pressure
91real :: dz ! altitude diff
92real :: signe ! orientation of lon axis for mountain torque computation
93
94include "planet.h"
95
96!===============================================================================
97! 1. Input parameters
98!===============================================================================
99
100pi = 2.*asin(1.)
101
102write(*,*) ""
103write(*,*) "You are working on the atmosphere of ",planet
104
105!===============================================================================
106! 1.1 Input file
107!===============================================================================
108
109write(*,*) ""
110write(*,*) " Program valid for Venus CAM or LMD output files"
111write(*,*) "Enter input file name:"
112
113read(*,'(a128)') infile
114write(*,*) ""
115
116! open input file
117
118ierr = NF_OPEN(infile,NF_NOWRITE,infid)
119if (ierr.ne.NF_NOERR) then
120   write(*,*) 'ERROR: Pb opening file ',trim(infile)
121   stop ""
122endif
123
124!===============================================================================
125! 1.2 Get grids in lon,lat,alt(pressure),time
126!===============================================================================
127
128call get_iddim(infid,lat_varid,latlength,lon_varid,lonlength,&
129                     alt_varid,altlength,time_varid,timelength,lmdflag )
130
131allocate(lon(lonlength))
132ierr=NF_GET_VAR_REAL(infid,lon_varid,lon)
133if (ierr.ne.NF_NOERR) stop "Error: Failed reading longitude"
134if(lon(1).gt.lon(2)) then
135  signe=-1.
136else
137  signe=1.
138endif
139
140allocate(lat(latlength))
141ierr=NF_GET_VAR_REAL(infid,lat_varid,lat)
142if (ierr.ne.NF_NOERR) stop "Error: Failed reading lat"
143
144allocate(latrad(latlength))
145latrad = lat*pi/180.
146
147! Lat, lon pressure intervals
148deltalat = abs(latrad(2)-latrad(1))
149deltalon = abs(lon(2)-lon(1))*pi/180.
150
151allocate(plev(altlength))
152ierr=NF_GET_VAR_REAL(infid,alt_varid,plev)
153if (ierr.ne.NF_NOERR) stop "Error: Failed reading pressure levels"
154if(.not.lmdflag) then
155! in CAM files, pressure in mbar and reversed...
156  call reverselev(altlength,plev)
157  plev=plev*100.  ! convertion to Pa
158endif
159
160allocate(time(timelength))
161ierr=NF_GET_VAR_REAL(infid,time_varid,time)
162if (ierr.ne.NF_NOERR) stop "Error: Failed reading time"
163
164! Time axis IN PLANET DAYS
165
166if(.not.lmdflag) then
167! in CAM files, time in Earth days...
168!   => seconds
169  time=time*86400.
170endif
171time=time/localday
172
173!===============================================================================
174! 1.3 Get output file name
175!===============================================================================
176write(*,*) ""
177!write(*,*) "Enter output file name"
178!read(*,*) outfile
179outfile=infile(1:len_trim(infile)-3)//"_LTRQ.nc"
180write(*,*) "Output file name is: "//trim(outfile)
181
182
183
184!===============================================================================
185! 2.1 Store needed fields
186!===============================================================================
187
188!===============================================================================
189! 2.1.1 Surface pressure and geopotential
190!===============================================================================
191allocate(ps(lonlength,latlength,timelength))
192allocate(phis3d(lonlength,latlength,timelength))
193allocate(phis(lonlength,latlength))
194
195text="PS"
196call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
197if (ierr1.ne.NF_NOERR) then
198  write(*,*) "  looking for psol instead... "
199  text="psol"
200  call get_var3d(infid,lonlength,latlength,timelength,text,ps,ierr1,ierr2)
201  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get psol ID"
202endif
203if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface pressure"
204if((.not.lmdflag).and.(planet.eq."Venus")) call reverse3d(lonlength,latlength,timelength,ps)
205
206text="PHIS"
207call get_var3d(infid,lonlength,latlength,timelength,text,phis3d,ierr1,ierr2)
208if (ierr1.ne.NF_NOERR) then
209  write(*,*) " Failed to get PHIS ID (3d), looking for phis (2d) instead... "
210  text="phis"
211  call get_var2d(infid,lonlength,latlength,text,phis,ierr1,ierr2)
212  if (ierr1.ne.NF_NOERR) stop "Error: Failed to get phis ID"
213  if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface geopotential"
214else
215  if (ierr2.ne.NF_NOERR) stop "Error: Failed reading surface geopotential"
216  phis(:,:)=phis3d(:,:,1)
217  if((.not.lmdflag).and.(planet.eq."Venus")) call reverse2d(lonlength,latlength,phis)
218endif
219
220if(lmdflag) then
221  p0=1.
222else
223  ierr=NF_INQ_VARID(infid,"P0",tmpvarid)
224  if (ierr.ne.NF_NOERR) then
225    write(*,*) "Failed to get P0 ID, used psref=",psref
226    p0=psref
227  else
228    ierr=NF_GET_VAR_REAL(infid,tmpvarid,p0)
229    if (ierr.ne.NF_NOERR) then
230      write(*,*) "Failed reading reference pressure, used psref=",psref
231      p0=psref
232    endif
233  endif
234endif
235
236!===============================================================================
237! 2.1.3 Winds
238!===============================================================================
239allocate(vitu(lonlength,latlength,altlength,timelength))
240
241! zonal wind vitu / U (in m/s)
242if(lmdflag) then
243  text="vitu"
244else
245  text="U"
246endif
247
248call get_var4d(infid,lonlength,latlength,altlength,timelength,text,vitu,ierr1,ierr2)
249if (ierr1.ne.NF_NOERR) stop "Error: Failed to get U ID"
250if (ierr2.ne.NF_NOERR) stop "Error: Failed reading zonal wind"
251
252if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,vitu)
253
254!===============================================================================
255! 2.1.5 Friction in Boundary layer
256!===============================================================================
257allocate(duvdf(lonlength,latlength,altlength,timelength))
258
259if(lmdflag) then
260  text="duvdf"
261else
262  text="DUVDF"
263endif
264call get_var4d(infid,lonlength,latlength,altlength,timelength,text,duvdf,ierr1,ierr2)
265if (ierr1.ne.NF_NOERR) then
266  write(*,*) "Failed to get duvdf ID"
267  flag_duvdf = .false.
268
269! IDEAL FRICTION ?
270write(*,*) ""
271write(*,*) " Is the friction at surface ideal ? 0 for no, 1 for yes"
272write(*,*) " duvdf = -u/(86400*30) in surface layer"
273write(*,*) " timescale hard-coded: 30 Edays"
274read(*,'(i1)') idlsurf
275!write(*,*) ""
276!write(*,*) " ASSUMED YES ! "
277!idlsurf=1
278write(*,*) ""
279
280if (idlsurf.eq.1) then
281  flag_duvdf = .true.
282  duvdf = 0.
283  ilev=1
284do ilon=1,lonlength
285 do ilat=1,latlength
286  do itim=1,timelength
287   duvdf(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) * (-1.)/(86400.*30.)
288  enddo
289 enddo
290enddo
291endif
292
293else
294  if (ierr2.ne.NF_NOERR) stop "Error: Failed reading duvdf"
295  if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,duvdf)
296  flag_duvdf = .true.
297endif
298
299!===============================================================================
300! 2.1.7 Dynamics (includes the mountain torque...)
301!===============================================================================
302allocate(dudyn(lonlength,latlength,altlength,timelength))
303
304if(lmdflag) then
305  text="dudyn"
306else
307  text="DUDYN"
308endif
309call get_var4d(infid,lonlength,latlength,altlength,timelength,text,dudyn,ierr1,ierr2)
310if (ierr1.ne.NF_NOERR) then
311  write(*,*) "Failed to get dudyn ID"
312  flag_dudyn = .false.
313else
314  if (ierr2.ne.NF_NOERR) stop "Error: Failed reading dudyn"
315  if(.not.lmdflag) call reverse4d(lonlength,latlength,altlength,timelength,dudyn)
316  flag_dudyn = .true.
317endif
318
319!===============================================================================
320! 2.2 Computations
321!===============================================================================
322
323!===============================================================================
324! 2.2.1 Pressure intervals
325!===============================================================================
326allocate(deltap(lonlength,latlength,altlength,timelength))
327allocate(press(lonlength,latlength,altlength,timelength))
328allocate(pp(altlength))
329
330if(lmdflag) then  ! LMD vs CAM
331  text="pres"
332  call get_var4d(infid,lonlength,latlength,altlength,timelength,text,press,ierr1,ierr2)
333  if (ierr1.ne.NF_NOERR) then
334    stop "Error: Failed to get pres ID"
335  else
336    if (ierr2.ne.NF_NOERR) stop "Error: Failed reading pres"
337  endif
338
339else  ! LMD vs CAM
340  ! hybrid coordinate aps/hyam
341  text="hyam"
342  ierr=NF_INQ_VARID(infid,text,tmpvarid)
343  if (ierr.ne.NF_NOERR) then
344    stop "Error: Could not find aps/hyam coordinates"
345  else
346    allocate(aps(altlength))
347    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
348    if (ierr.ne.NF_NOERR) then
349      stop "Error: Failed reading aps/hyam"
350    endif
351    call reverselev(altlength,aps)
352  endif
353
354  ! hybrid coordinate hybm
355  text="hybm"
356  ierr=NF_INQ_VARID(infid,text,tmpvarid)
357  if (ierr.ne.NF_NOERR) then
358    stop "Error: Failed to get bps/hybm ID"
359  else
360    allocate(bps(altlength))
361    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
362    if (ierr.ne.NF_NOERR) then
363      stop "Error: Failed reading bps/hybm"
364    endif
365    call reverselev(altlength,bps)
366  endif
367
368  aps=aps*p0
369
370  do itim=1,timelength
371    do ilev=1,altlength
372      do ilat=1,latlength
373        do ilon=1,lonlength
374          press(ilon,ilat,ilev,itim)=aps(ilev)+bps(ilev)*ps(ilon,ilat,itim)
375        enddo
376      enddo
377    enddo
378  enddo
379
380endif  ! LMD vs CAM
381
382do itim=1,timelength
383
384do ilon=1,lonlength
385 do ilat=1,latlength
386  tmpp=ps(ilon,ilat,itim)
387  pp=press(ilon,ilat,:,itim)
388  deltap(ilon,ilat,1,itim)= tmpp - exp((log(pp(1))+log(pp(2)))/2.) ! initialization, rectified with ps below
389do ilev=2,altlength-1
390  deltap(ilon,ilat,ilev,itim)= exp((log(pp(ilev-1))+log(pp(ilev)))/2.)-&
391                     exp((log(pp(ilev))+log(pp(ilev+1)))/2.)
392enddo
393deltap(ilon,ilat,altlength,itim)=exp((log(pp(altlength-1))+log(pp(altlength)))/2.)
394 enddo
395enddo
396
397enddo ! timelength
398
399!===============================================================================
400! 2.2.2 Mass in cells
401!===============================================================================
402allocate(grav(lonlength,latlength,altlength,timelength))
403allocate(dmass(lonlength,latlength,altlength,timelength))
404
405do itim=1,timelength
406
407do ilon=1,lonlength
408 do ilat=1,latlength
409  do ilev=1,altlength
410    if (deltap(ilon,ilat,ilev,itim).ne.miss_val) then
411! Need to be consistent with GCM computations
412      grav(ilon,ilat,ilev,itim) = g0
413     dmass(ilon,ilat,ilev,itim) = &
414               a0*cos(latrad(ilat))*deltalon &
415             * a0*deltalat &
416             * deltap(ilon,ilat,ilev,itim)/grav(ilon,ilat,ilev,itim)
417    else
418      grav(ilon,ilat,ilev,itim) = miss_val
419     dmass(ilon,ilat,ilev,itim) = miss_val
420    endif
421  enddo
422 enddo
423enddo
424
425enddo ! timelength
426
427!===============================================================================
428! 2.2.3 Specific angular momentum
429!===============================================================================
430allocate(osam(lonlength,latlength,altlength,timelength))
431allocate(rsam(lonlength,latlength,altlength,timelength))
432
433do itim=1,timelength
434
435do ilon=1,lonlength
436 do ilat=1,latlength
437  do ilev=1,altlength
438      osam(ilon,ilat,ilev,itim) = &
439         a0*a0* cos(latrad(ilat))*cos(latrad(ilat)) * omega
440      rsam(ilon,ilat,ilev,itim) = vitu(ilon,ilat,ilev,itim) &
441       * a0*cos(latrad(ilat))
442  enddo
443 enddo
444enddo
445
446enddo ! timelength
447
448!===============================================================================
449! 2.2.4 Total angular momentum
450!===============================================================================
451allocate(oaam(timelength))
452allocate(raam(timelength))
453
454do itim=1,timelength
455
456oaam(itim) = 0.
457raam(itim) = 0.
458do ilon=1,lonlength
459 do ilat=1,latlength
460  do ilev=1,altlength
461      oaam(itim) = oaam(itim) &
462       + osam(ilon,ilat,ilev,itim)/ hadday * dmass(ilon,ilat,ilev,itim)
463      raam(itim) = raam(itim) &
464       + rsam(ilon,ilat,ilev,itim)/ hadday * dmass(ilon,ilat,ilev,itim)
465  enddo
466 enddo
467enddo
468if (oaam(itim).eq.0.) then
469  oaam(itim) = miss_val
470  raam(itim) = miss_val
471endif
472
473enddo ! timelength
474
475!===============================================================================
476! 2.2.5 Mountain, friction, dynamics and GW torques
477!===============================================================================
478allocate(tmou(latlength,timelength))
479if (flag_dudyn)  allocate(tdyn(latlength,timelength))
480if (flag_duvdf)  allocate(tbls(latlength,timelength))
481
482do itim=1,timelength
483
484tmou(:,itim) = 0.
485 do ilat=2,latlength-1
486do ilon=1,lonlength
487      if (ilon.eq.lonlength) then
488         dz = (phis(   1,ilat)     -phis(ilon,ilat))/g0
489       tmpp = (ps(     1,ilat,itim)  +ps(ilon,ilat,itim))/2.
490      else
491         dz = (phis(ilon+1,ilat)   -phis(ilon,ilat))/g0
492       tmpp = (ps(ilon+1,ilat,itim)  +ps(ilon,ilat,itim))/2.
493      endif
494      tmou(ilat,itim) = tmou(ilat,itim) + a0*a0* deltalat*cos(latrad(ilat)) &
495       * signe*dz*tmpp &
496       / hadley
497enddo
498 enddo
499
500if (flag_dudyn) then
501tdyn(:,itim) = 0.
502 do ilat=1,latlength
503do ilon=1,lonlength
504  do ilev=1,altlength
505      tdyn(ilat,itim) = tdyn(ilat,itim) + dudyn(ilon,ilat,ilev,itim) &
506       * a0*cos(latrad(ilat)) * dmass(ilon,ilat,ilev,itim) &
507       / hadley
508  enddo
509enddo
510  if (tdyn(ilat,itim).eq.0.) tdyn(ilat,itim) = miss_val
511 enddo
512endif
513
514if (flag_duvdf) then
515tbls(:,itim) = 0.
516 do ilat=1,latlength
517do ilon=1,lonlength
518  do ilev=1,altlength
519      tbls(ilat,itim) = tbls(ilat,itim) + duvdf(ilon,ilat,ilev,itim) &
520       * a0*cos(latrad(ilat)) * dmass(ilon,ilat,ilev,itim) &
521       / hadley
522  enddo
523enddo
524  if (tbls(ilat,itim).eq.0.) tbls(ilat,itim) = miss_val
525 enddo
526endif
527
528enddo ! timelength
529
530print*,"End of computations"
531
532!===============================================================================
533! 3. Create output file
534!===============================================================================
535
536! Create output file
537ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
538if (ierr.ne.NF_NOERR) then
539  write(*,*)"Error: could not create file ",outfile
540  stop
541endif
542
543!===============================================================================
544! 3.1. Define and write dimensions
545!===============================================================================
546
547call write_dim(outfid,lonlength,latlength,altlength,timelength, &
548    lon,lat,plev,time,lon_dimid,lat_dimid,alt_dimid,time_dimid)
549
550!===============================================================================
551! 3.2. Define and write variables
552!===============================================================================
553
554! 1D Variables
555
556datashape1d=time_dimid
557 
558call write_var1d(outfid,datashape1d,timelength,&
559                "oaam      ", "Total rotation ang  ","E25kgm2s-1",miss_val,&
560                 oaam )
561
562call write_var1d(outfid,datashape1d,timelength,&
563                "raam      ", "Total wind ang      ","E25kgm2s-1",miss_val,&
564                 raam )
565
566! 2D variables
567
568datashape2d(1)=lat_dimid
569datashape2d(2)=time_dimid
570
571call write_var2d(outfid,datashape2d,latlength,timelength,&
572                "tmou      ", "Mountain torque     ","E18kgm2s-2",miss_val,&
573                 tmou )
574
575if (flag_dudyn) then
576call write_var2d(outfid,datashape2d,latlength,timelength,&
577                "tdyn      ", "Dynamics torque     ","E18kgm2s-2",miss_val,&
578                 tdyn )
579endif
580
581if (flag_duvdf) then
582call write_var2d(outfid,datashape2d,latlength,timelength,&
583                "tbls      ", "Friction torque     ","E18kgm2s-2",miss_val,&
584                 tbls )
585endif
586
587datashape2d(1)=lon_dimid
588datashape2d(2)=lat_dimid
589
590call write_var2d(outfid,datashape2d,lonlength,latlength,&
591                "phis      ", "Surface geopot      ","m2 s-2    ",miss_val,&
592                 phis )
593
594! 3D variables
595
596datashape3d(1)=lon_dimid
597datashape3d(2)=lat_dimid
598datashape3d(3)=time_dimid
599
600call write_var3d(outfid,datashape3d,lonlength,latlength,timelength,&
601                "ps        ", "Surface pressure    ","Pa        ",miss_val,&
602                 ps )
603
604datashape3d(1)=lat_dimid
605datashape3d(2)=alt_dimid
606datashape3d(3)=time_dimid
607
608!!!! Close output file
609ierr=NF_CLOSE(outfid)
610if (ierr.ne.NF_NOERR) then
611  write(*,*) 'Error, failed to close output file ',outfile
612endif
613
614
615end program
616
617subroutine reverse4d(nlon,nlat,np,nt,newf)
618!==============================================================================
619! Purpose:
620! reverse lat, lon and vertical axis on a 3D+time field
621!==============================================================================
622implicit none
623!==============================================================================
624! Arguments:
625!==============================================================================
626integer,intent(in) :: nlon ! longitude size
627integer,intent(in) :: nlat ! latitude size
628integer,intent(in) :: np   ! vertical size
629integer,intent(in) :: nt   ! time size
630real,intent(inout),dimension(nlon,nlat,np,nt) :: newf ! field to be reversed
631
632!==============================================================================
633! Local variables:
634!==============================================================================
635integer :: i,j,k,l
636real,dimension(nlon,nlat,np,nt) :: tmp,tmp2 ! reversed fields
637include "planet.h"
638
639! Vertical axis:
640    do k=1,np
641          tmp(:,:,k,:)=newf(:,:,np+1-k,:)
642    enddo
643! horizontal dimensions
644if(planet.eq."Venus") then
645  do l=1,nt
646      do i=1,nlat
647        do j=1,nlon
648          tmp2(j,i,:,l)=tmp(nlon+1-j,nlat+1-i,:,l)
649        enddo
650      enddo
651  enddo
652else
653  tmp2=tmp
654endif
655
656newf=tmp2
657
658end subroutine reverse4d
659
660!===============================================================================
661
662subroutine reverse2d(nlon,nlat,newf)
663!==============================================================================
664! Purpose:
665! reverse lat and lon on a 2D field
666!==============================================================================
667implicit none
668!==============================================================================
669! Arguments:
670!==============================================================================
671integer,intent(in) :: nlon ! longitude size
672integer,intent(in) :: nlat ! latitude size
673real,intent(inout),dimension(nlon,nlat) :: newf ! field to be reversed
674
675!==============================================================================
676! Local variables:
677!==============================================================================
678integer :: i,j
679real,dimension(nlon,nlat) :: tmp ! reversed field
680
681      do i=1,nlat
682        do j=1,nlon
683          tmp(j,i)=newf(nlon+1-j,nlat+1-i)
684        enddo
685      enddo
686newf=tmp
687
688end subroutine reverse2d
689
690!===============================================================================
691
692subroutine reverse3d(nlon,nlat,nt,newf)
693!==============================================================================
694! Purpose:
695! reverse lat and lon on a 2D+time field
696!==============================================================================
697implicit none
698!==============================================================================
699! Arguments:
700!==============================================================================
701integer,intent(in) :: nlon ! longitude size
702integer,intent(in) :: nlat ! latitude size
703integer,intent(in) :: nt ! time size
704real,intent(inout),dimension(nlon,nlat,nt) :: newf ! field to be reversed
705
706!==============================================================================
707! Local variables:
708!==============================================================================
709integer :: i,j,l
710real,dimension(nlon,nlat,nt) :: tmp ! reversed field
711
712  do l=1,nt
713      do i=1,nlat
714        do j=1,nlon
715          tmp(j,i,l)=newf(nlon+1-j,nlat+1-i,l)
716        enddo
717      enddo
718  enddo
719newf=tmp
720
721end subroutine reverse3d
722
723!===============================================================================
724
725subroutine reverselev(np,newf)
726!==============================================================================
727! Purpose:
728! reverse vertical pressure axis
729!==============================================================================
730implicit none
731!==============================================================================
732! Arguments:
733!==============================================================================
734integer,intent(in) :: np   ! vertical size
735real,intent(inout),dimension(np) :: newf ! field to be reversed
736
737!==============================================================================
738! Local variables:
739!==============================================================================
740integer :: k
741real,dimension(np) :: tmp ! reversed field
742
743    do k=1,np
744          tmp(k)=newf(np+1-k)
745    enddo
746newf=tmp
747
748end subroutine reverselev
Note: See TracBrowser for help on using the repository browser.