Changeset 414

Timestamp:

06/10/16 20:36:37 (8 years ago)

Author:

ymipsl

Message:

Etat0 of dcmip2016 cyclone use the given initialisation subroutine instaed of dcmip2012 initialisation

YM

Location:

codes/icosagcm/trunk/src

Files:

• dcmip2016_cyclone.f90 (added)
• etat0_dcmip2016_cyclone.f90 (modified) (4 diffs)

codes/icosagcm/trunk/src/etat0_dcmip2016_cyclone.f90

@@ -3 +3 @@
   IMPLICIT NONE
   PRIVATE
-  REAL(rstd),PARAMETER :: zt=15000
-  REAL(rstd),PARAMETER :: q0=0.021
-  REAL(rstd),PARAMETER :: qt=1e-11
-  REAL(rstd),PARAMETER :: T0=302.15
-  REAL(rstd),PARAMETER :: Ts=302.15
-  REAL(rstd),PARAMETER :: zq1=3000
-  REAL(rstd),PARAMETER :: zq2=8000
-  REAL(rstd),PARAMETER :: Gamma=0.007
-  REAL(rstd),PARAMETER :: pb=101500
-  REAL(rstd),PARAMETER :: Deltap=1115
-  REAL(rstd),PARAMETER :: rp=282000
-  REAL(rstd),PARAMETER :: zp=7000
-  REAL(rstd),PARAMETER :: epsilon=1e-25
-  REAL(rstd),PARAMETER :: Rd=287
-  
-  REAL(rstd), PARAMETER :: lonc=Pi, latc=Pi/18, &
-       Tv0=T0*(1+0.608*q0), Tvt=Tv0-Gamma*zt
+
 
   PUBLIC getin_etat0, compute_etat0
@@ -28 +12 @@
     USE mpipara, ONLY : is_mpi_root
 
-    IF(nqtot<1) THEN
+    IF(nqtot<5) THEN
        IF (is_mpi_root)  THEN
-          PRINT *, "nqtot must be at least 1 for test case DCMIP5"
+          PRINT *, "nqtot must be at least 5 for test case DCMIP2016 tropical cyclone"
        END IF
        STOP
@@ -38 +22 @@
 
   SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, Temp, ulon, ulat, q)
-    USE disvert_mod
-    USE omp_para
+  USE disvert_mod
+  USE omp_para
+  USE dcmip2016_cyclone_mod, ONLY : tropical_cyclone_test
+  USE terminator, ONLY: initial_value_Terminator
+  IMPLICIT NONE
     INTEGER, INTENT(IN)    :: ngrid
     REAL(rstd),INTENT(IN)  :: lon(ngrid)
@@ -50 +37 @@
     REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
 
+    INTEGER :: zcoords
+    REAL :: p,z
+    REAL :: rho, thetav
     INTEGER :: l,ij
-    REAL(rstd) :: r, d, d1,d2, zz, Tave, Tp, vt, aa, bb
 
-    phis(:)=0.
-
+   
     DO ij=1,ngrid
-       r=radius*arc(lonc,latc, lon(ij),lat(ij))
-       ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5)
-    END DO
-
-    DO l=ll_begin,ll_end
-       aa=.5*(ap(l)+ap(l+1))
-       bb=.5*(bp(l)+bp(l+1))
-       DO ij=1,ngrid
-          r=radius*arc(lonc,latc, lon(ij),lat(ij))
-          CALL compute_z(aa+bb*ps(ij),ps(ij),r,zz)
-          IF (zz<=zt) THEN
-             q(ij,l,1)=q0*exp(-zz/zq1)*exp(-(zz/zq2)**2)
-             Tave=Tv0-Gamma*zz
-             Tp=(Tv0-Gamma*zz) * ( ( 1 + 2*Rd*(Tv0-Gamma*zz)*zz /                          &
-                  (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(zz/zp)**2))))**(-1) - 1) 
-             vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
-                  / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
-                  - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
-          ELSE
-             q(ij,l,1)=qt
-             Tave=Tvt
-             Tp=0
-             vt=0
-          ENDIF
-          Temp(ij,l)=Tave+Tp 
-
-          d1=sin(latc)*cos(lat(ij))-cos(latc)*sin(lat(ij))*cos(lon(ij)-lonc)
-          d2=cos(latc)*sin(lon(ij)-lonc)
-          d=MAX(1e-25,sqrt(d1**2+d2**2))
-          ulon(ij,l)=vt*d1/d
-          ulat(ij,l)=vt*d2/d
+       z=0.
+       zcoords=1                 
+       CALL tropical_cyclone_test(lon(ij),lat(ij),p,z,zcoords,ulon(ij,1),ulat(ij,1), &
+                                 temp(ij,1),thetav,phis(ij),ps(ij),rho,q(ij,1,1))
+       
+       zcoords=0
+       DO l=ll_begin,ll_end
+         p=0.5*(ap(l)+ap(l+1) + (bp(l)+bp(l+1)) * ps(ij))
+         CALL tropical_cyclone_test(lon(ij),lat(ij),p,z,zcoords,ulon(ij,l),ulat(ij,l), &
+                                   temp(ij,l),thetav,phis(ij),ps(ij),rho,q(ij,l,1))
+                                   
+          IF (physics_thermo==thermo_fake_moist) temp(ij,l)=Temp(ij,l)*(1+0.608*q(ij,l,1))
+          q(ij,l,2)=0.
+          q(ij,l,3)=0.
+          CALL initial_value_Terminator(lat(ij),lon(ij),q(ij,l,4),q(ij,l,5))
        END DO
     ENDDO
+
+
   END SUBROUTINE compute_etat0
   
-  SUBROUTINE compute_z(pmodel,ps,r,z)
-  USE icosa
-  IMPLICIT NONE
-     REAL(rstd),PARAMETER   :: epsilon0=2e-13
-     REAL(rstd),INTENT(IN)  :: pmodel
-     REAL(rstd),INTENT(IN)  :: ps
-     REAL(rstd),INTENT(IN)  :: r
-     REAL(rstd),INTENT(OUT) :: z
-     
-     REAL(rstd) :: pt, pave, pp, dpdz
-     REAL(rstd) :: znp1
-     REAL(rstd) :: epsilon
-     REAL(rstd) :: p
-     INTEGER  :: n     
 
-     pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma))
-
-     IF (pmodel>pt) THEN
-       z=Tv0/Gamma*(1-(pmodel/ps)**(Rd*Gamma/g))
-     ELSE
-       z=zt+Rd*Tvt/g*log(Pt/pmodel)
-     ENDIF
-     
-     IF (z>zt .OR. r>1000000) return
-     
-     epsilon=1
-     n=0
-     
-     DO WHILE(epsilon>epsilon0 .AND. n<20)
-  
-       IF (z<zt) THEN
-         pave=pb*((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))
-         pp= -Deltap * exp(-(r/rp)**1.5-(z/zp)**2) *  ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma))
-       ELSE
-         pave=pt*exp(g*(zt-z)/(Rd*Tvt))
-         pp=0
-       ENDIF
-       p=pave+pp
-     
-       dpdz =  2*Deltap*z/zp**2 * exp(-(r/rp)**1.5) * exp(-(z/zp)**2) * ((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))  &
-             - g/(Rd*Tv0) * (preff - Deltap*exp(-(r/rp)**1.5) * exp(-(z/zp)**2)) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)-1)
-           
-       znp1 = z - ( pmodel - p ) / (-dpdz)
-       
-       epsilon=ABS( (znp1-z)/znp1 )
-       z=znp1
-       n=n+1
-     ENDDO  
-   
-   END SUBROUTINE compute_z
 
 END MODULE etat0_dcmip2016_cyclone_mod