Changeset 1442 for trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/newstart.F

Timestamp:

Jun 4, 2015, 4:23:32 PM (10 years ago)

Author:

slebonnois

Message:

SL: update of the Venus GCM, + corrections on routines used for newstart/start2archive for Titan and Venus, + some modifications on tools

File:

• trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/newstart.F (modified) (10 diffs)

trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/newstart.F

@@ -24 +24 @@
       use exner_hyb_m, only: exner_hyb
       use exner_milieu_m, only: exner_milieu
+      USE comconst_mod
+      USE comvert_mod, ONLY: ap,bp,presnivs,pa,preff,nivsigs,nivsig,
+     .                       aps,bps,scaleheight,pseudoalt,
+     .                       disvert_type,pressure_exner
 
       implicit none
@@ -29 +33 @@
 #include "dimensions.h"
 #include "paramet.h"
-#include "comconst.h"
 #include "comdissnew.h"
-#include "comvert.h"
 #include "comgeom2.h"
 #include "logic.h"
@@ -140 +142 @@
       integer, dimension(4) :: start,counter
       REAL phisinverse(iip1,jjp1)  ! geopotentiel au sol avant inversion
-      logical topoflag,albedoflag,razvitu,razvitv 
+      logical topoflag,notopo,albedoflag,razvitu,razvitv,uini
+      logical razTS,raztemp 
+      real, dimension(:), allocatable :: tvira,dzst,zkm
       real    albedo
       
@@ -880 +884 @@
        topoflag = . FALSE .
        CALL getin('topoflag',topoflag)
+! notopo = T: we go back to flat surface
+       notopo = .FALSE.
+       CALL getin('notopo',notopo)
 
         print*,zmeaold(2,1:10)
@@ -906 +913 @@
      .            0.0,jjm,rlonu,rlatv,.true.)
 
+       ELSE IF ( notopo ) THEN
+          print*,'Flattening the topography'
+          phis=0.
+          zmea=0.
+          zstd=0.
+          zsig=0.
+          zgam=0.
+          zthe=0.
+          zpic=0.
+          zval=0.
        ELSE
           print*,'Using existing topography'
@@ -950 +967 @@
 
 c Temperature de surface
-      call interp_horiz (tsurfold,tsurfS,imold,jmold,iim,jjm,1,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (1,iip1,jjp1,ngridmx,tsurfS,tsurf)
-c     write(44,*) 'tsurf', tsurf
+! razTS need to be in the specific run.def for newstart
+      razTS = . FALSE .
+      CALL getin('razTS',razTS)
+
+      if (razTS) then
+        tsurf(:) = 735.
+      else  
+       call interp_horiz (tsurfold,tsurfS,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+       call gr_dyn_fi (1,iip1,jjp1,ngridmx,tsurfS,tsurf)
+c      write(44,*) 'tsurf', tsurf
+      endif
 
 c Temperature du sous-sol
-      call interp_horiz(tsoilold,tsoilS,
+      if (razTS) then
+        tsoil(:,:)=735.
+      else  
+       call interp_horiz(tsoilold,tsoilS,
      &                  imold,jmold,iim,jjm,nsoilmx,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (nsoilmx,iip1,jjp1,ngridmx,tsoilS,tsoil)
-c     write(45,*) 'tsoil',tsoil
+       call gr_dyn_fi (nsoilmx,iip1,jjp1,ngridmx,tsoilS,tsoil)
+c      write(45,*) 'tsoil',tsoil
+      endif
 
 ! CHANGING ALBEDO: may be done through run.def
@@ -1050 +1079 @@
       CALL pression(ip1jmp1, ap, bp, ps, p3d)
          if (disvert_type==1) then
-           CALL exner_hyb(  ip1jmp1, ps, p3d, pks, pk, pkf )
+           CALL exner_hyb(  ip1jmp1, ps, p3d,pks, pk, pkf )
          else ! we assume that we are in the disvert_type==2 case
            CALL exner_milieu( ip1jmp1, ps, p3d, pks, pk, pkf )
@@ -1067 +1096 @@
 c     enddo
 
+! raztemp need to be in the specific run.def for newstart
+      raztemp = . FALSE .
+      CALL getin('raztemp',raztemp)
+
+! Reinitialisation of temperature to VIRA profile lisse
+      if (raztemp) then
+
+        allocate(tvira(0:lmold),dzst(0:lmold),zkm(0:lmold))
+        print*,"Venus = temperature initiale imposee = VIRA lisse "
+        dzst(0) = 0.0
+        dzst(1) = -log(p3d(1,1,2)/preff)*r/g
+        do l=2,lmold
+           dzst(l)=-log(p3d(1,1,l+1)/p3d(1,1,l))*r/g
+        enddo
+        tvira(0) = 735.
+        zkm(0) = 0.0
+        do l=1,lmold
+          zkm(l) = zkm(l-1)+tvira(l-1)*dzst(l)/1000. ! approx avec T(l-1)
+          if(zkm(l).lt.60.) then
+            tvira(l)=735.-7.95*zkm(l)
+          else
+            tvira(l)=AMAX1(258.-3.*(zkm(l)-60.),168.)
+          endif
+          zkm(l) = zkm(l-1)+(tvira(l-1)+tvira(l))/2.*dzst(l)/1000.
+        enddo
+        do l=1,lmold
+         do j=1,jmold+1
+          do i=1,imold+1
+            told(i,j,l)=tvira(l)
+          enddo
+         enddo
+        enddo
+      endif  ! end raztemp
+
       write (*,*) 'told ', told (1,jmold+1,1)  ! INFO
       call interp_vert
@@ -1087 +1150 @@
       teta(iip1,:,:) =  teta(1,:,:)
 
+! RESETING U TO uini: may be done through run.def
+       uini = .FALSE.
+       CALL getin('uini',uini)
 ! RESETING U TO 0: may be done through run.def
        razvitu = . FALSE .
@@ -1110 +1176 @@
       call scal_wind(us,vs,unat,vnat)
 ! Reseting u=0
-      if (razvitu) then
-           unat(:,:,:) = 0.
+      if ((razvitu).and..not.(uini)) then
+         unat(:,:,:) = 0.
+      endif
+! Reseting u=uini
+      if ((uini).and..not.(razvitu)) then
+         do j=1,jjp1
+           do l=1,llm
+             if (p3d(1,j,l).gt.3e3) then
+               unat(:,j,l) = -110./8.03*log(p3d(:,j,l)/9.2e6)
+             else
+               unat(:,j,l) = 110./6.62*(log(p3d(:,j,l)/9.2e6)+14.65)
+             endif
+             if (abs(rlatS(1,j)).gt.50.) then
+              unat(:,j,l)=unat(:,j,l)*(90.-abs(rlatS(:,j)))/40.
+             endif
+           enddo
+         enddo
+      endif
+! incompatible options
+      if ((uini).and.(razvitu)) then
+         print*,"You have to choose between razvitu and uini..."
+         stop
       endif
       write (*,*) 'unat ', unat (1,2,1)    ! INFO
+
       do l=1,llm
         do j = 1, jjp1