Changeset 5246 for LMDZ6/trunk/libf/dyn3d/sw_case_williamson91_6.f90
- Timestamp:
- Oct 21, 2024, 2:58:45 PM (32 hours ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3d/sw_case_williamson91_6.f90
r5245 r5246 2 2 ! $Id $ 3 3 ! 4 4 SUBROUTINE sw_case_williamson91_6(vcov,ucov,teta,masse,ps) 5 5 6 c======================================================================= 7 c 8 c Author: Thomas Dubos original: 26/01/2010 9 c ------- 10 c 11 c Subject: 12 c ------ 13 c Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz 14 c 15 c Method: 16 c -------- 17 c 18 c Interface: 19 c ---------- 20 c 21 c Input: 22 c ------ 23 c 24 c Output: 25 c ------- 26 c 27 c======================================================================= 28 USE comconst_mod, ONLY: cpp, omeg, rad 29 USE comvert_mod, ONLY: ap, bp, preff 30 31 IMPLICIT NONE 32 c----------------------------------------------------------------------- 33 c Declararations: 34 c --------------- 6 !======================================================================= 7 ! 8 ! Author: Thomas Dubos original: 26/01/2010 9 ! ------- 10 ! 11 ! Subject: 12 ! ------ 13 ! Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz 14 ! 15 ! Method: 16 ! -------- 17 ! 18 ! Interface: 19 ! ---------- 20 ! 21 ! Input: 22 ! ------ 23 ! 24 ! Output: 25 ! ------- 26 ! 27 !======================================================================= 28 USE comconst_mod, ONLY: cpp, omeg, rad 29 USE comvert_mod, ONLY: ap, bp, preff 35 30 36 include "dimensions.h"37 include "paramet.h"38 include "comgeom.h"39 include "iniprint.h"31 IMPLICIT NONE 32 !----------------------------------------------------------------------- 33 ! Declararations: 34 ! --------------- 40 35 41 c Arguments: 42 c ---------- 36 include "dimensions.h" 37 include "paramet.h" 38 include "comgeom.h" 39 include "iniprint.h" 43 40 44 c variables dynamiques 45 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 46 REAL teta(ip1jmp1,llm) ! temperature potentielle 47 REAL ps(ip1jmp1) ! pression au sol 48 REAL masse(ip1jmp1,llm) ! masse d'air 49 REAL phis(ip1jmp1) ! geopotentiel au sol 41 ! Arguments: 42 ! ---------- 50 43 51 c Local: 52 c ------ 44 ! variables dynamiques 45 REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 46 REAL :: teta(ip1jmp1,llm) ! temperature potentielle 47 REAL :: ps(ip1jmp1) ! pression au sol 48 REAL :: masse(ip1jmp1,llm) ! masse d'air 49 REAL :: phis(ip1jmp1) ! geopotentiel au sol 53 50 54 REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches 55 REAL pks(ip1jmp1) ! exner au sol 56 REAL pk(ip1jmp1,llm) ! exner au milieu des couches 57 REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches 58 REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm) 51 ! Local: 52 ! ------ 59 53 60 REAL :: sinth,costh,costh2, Ath,Bth,Cth, lon,dps 61 INTEGER i,j,ij 54 REAL :: p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches 55 REAL :: pks(ip1jmp1) ! exner au sol 56 REAL :: pk(ip1jmp1,llm) ! exner au milieu des couches 57 REAL :: pkf(ip1jmp1,llm) ! exner filt.au milieu des couches 58 REAL :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm) 62 59 63 REAL, PARAMETER :: rho=1 ! masse volumique de l'air (arbitraire) 64 REAL, PARAMETER :: K = 7.848e-6 ! K = \omega 65 REAL, PARAMETER :: gh0 = 9.80616 * 8e3 66 INTEGER, PARAMETER :: R0=4, R1=R0+1, R2=R0+2 ! mode 4 67 c NB : rad = 6371220 dans W91 (6371229 dans LMDZ) 68 c omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ) 69 70 IF(0==0) THEN 71 c Williamson et al. (1991) : onde de Rossby-Haurwitz 72 teta = preff/rho/cpp 73 c geopotentiel (pression de surface) 74 do j=1,jjp1 75 costh2 = cos(rlatu(j))**2 76 Ath = (R0+1)*(costh2**2) + (2*R0*R0-R0-2)*costh2 - 2*R0*R0 77 Ath = .25*(K**2)*(costh2**(R0-1))*Ath 78 Ath = .5*K*(2*omeg+K)*costh2 + Ath 79 Bth = (R1*R1+1)-R1*R1*costh2 80 Bth = 2*(omeg+K)*K/(R1*R2) * (costh2**(R0/2))*Bth 81 Cth = R1*costh2 - R2 82 Cth = .25*K*K*(costh2**R0)*Cth 83 do i=1,iip1 84 ij=(j-1)*iip1+i 85 lon = rlonv(i) 86 dps = Ath + Bth*cos(R0*lon) + Cth*cos(2*R0*lon) 87 ps(ij) = rho*(gh0 + (rad**2)*dps) 88 enddo 89 enddo 90 write(lunout,*) 'W91 ps', MAXVAL(ps), MINVAL(ps) 91 c vitesse zonale ucov 92 do j=1,jjp1 93 costh = cos(rlatu(j)) 94 costh2 = costh**2 95 Ath = rad*K*costh 96 Bth = R0*(1-costh2)-costh2 97 Bth = rad*K*Bth*(costh**(R0-1)) 98 do i=1,iip1 99 ij=(j-1)*iip1+i 100 lon = rlonu(i) 101 ucov(ij,1) = (Ath + Bth*cos(R0*lon)) 102 enddo 103 enddo 104 write(lunout,*) 'W91 u', MAXVAL(ucov(:,1)), MINVAL(ucov(:,1)) 105 ucov(:,1)=ucov(:,1)*cu 106 c vitesse meridienne vcov 107 do j=1,jjm 108 sinth = sin(rlatv(j)) 109 costh = cos(rlatv(j)) 110 Ath = -rad*K*R0*sinth*(costh**(R0-1)) 111 do i=1,iip1 112 ij=(j-1)*iip1+i 113 lon = rlonv(i) 114 vcov(ij,1) = Ath*sin(R0*lon) 115 enddo 116 enddo 117 write(lunout,*) 'W91 v', MAXVAL(vcov(:,1)), MINVAL(vcov(:,1)) 118 vcov(:,1)=vcov(:,1)*cv 119 120 c ucov=0 121 c vcov=0 122 ELSE 123 c test non-tournant, onde se propageant en latitude 124 do j=1,jjp1 125 do i=1,iip1 126 ij=(j-1)*iip1+i 127 ps(ij) = 1e5*(1 + .1*exp(-100*(1+sin(rlatu(j)))**2) ) 128 enddo 129 enddo 130 131 c rho = preff/(cpp*teta) 132 teta = .01*preff/cpp ! rho = 100 ; phi = ps/rho = 1e3 ; c=30 m/s = 2600 km/j = 23 degres / j 133 ucov=0. 134 vcov=0. 135 END IF 136 137 CALL pression ( ip1jmp1, ap, bp, ps, p ) 138 CALL massdair(p,masse) 60 REAL :: sinth,costh,costh2, Ath,Bth,Cth, lon,dps 61 INTEGER :: i,j,ij 139 62 140 END 141 c----------------------------------------------------------------------- 63 REAL, PARAMETER :: rho=1 ! masse volumique de l'air (arbitraire) 64 REAL, PARAMETER :: K = 7.848e-6 ! K = \omega 65 REAL, PARAMETER :: gh0 = 9.80616 * 8e3 66 INTEGER, PARAMETER :: R0=4, R1=R0+1, R2=R0+2 ! mode 4 67 ! NB : rad = 6371220 dans W91 (6371229 dans LMDZ) 68 ! omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ) 69 70 IF(0==0) THEN 71 ! Williamson et al. (1991) : onde de Rossby-Haurwitz 72 teta = preff/rho/cpp 73 ! geopotentiel (pression de surface) 74 do j=1,jjp1 75 costh2 = cos(rlatu(j))**2 76 Ath = (R0+1)*(costh2**2) + (2*R0*R0-R0-2)*costh2 - 2*R0*R0 77 Ath = .25*(K**2)*(costh2**(R0-1))*Ath 78 Ath = .5*K*(2*omeg+K)*costh2 + Ath 79 Bth = (R1*R1+1)-R1*R1*costh2 80 Bth = 2*(omeg+K)*K/(R1*R2) * (costh2**(R0/2))*Bth 81 Cth = R1*costh2 - R2 82 Cth = .25*K*K*(costh2**R0)*Cth 83 do i=1,iip1 84 ij=(j-1)*iip1+i 85 lon = rlonv(i) 86 dps = Ath + Bth*cos(R0*lon) + Cth*cos(2*R0*lon) 87 ps(ij) = rho*(gh0 + (rad**2)*dps) 88 enddo 89 enddo 90 write(lunout,*) 'W91 ps', MAXVAL(ps), MINVAL(ps) 91 ! vitesse zonale ucov 92 do j=1,jjp1 93 costh = cos(rlatu(j)) 94 costh2 = costh**2 95 Ath = rad*K*costh 96 Bth = R0*(1-costh2)-costh2 97 Bth = rad*K*Bth*(costh**(R0-1)) 98 do i=1,iip1 99 ij=(j-1)*iip1+i 100 lon = rlonu(i) 101 ucov(ij,1) = (Ath + Bth*cos(R0*lon)) 102 enddo 103 enddo 104 write(lunout,*) 'W91 u', MAXVAL(ucov(:,1)), MINVAL(ucov(:,1)) 105 ucov(:,1)=ucov(:,1)*cu 106 ! vitesse meridienne vcov 107 do j=1,jjm 108 sinth = sin(rlatv(j)) 109 costh = cos(rlatv(j)) 110 Ath = -rad*K*R0*sinth*(costh**(R0-1)) 111 do i=1,iip1 112 ij=(j-1)*iip1+i 113 lon = rlonv(i) 114 vcov(ij,1) = Ath*sin(R0*lon) 115 enddo 116 enddo 117 write(lunout,*) 'W91 v', MAXVAL(vcov(:,1)), MINVAL(vcov(:,1)) 118 vcov(:,1)=vcov(:,1)*cv 119 120 ! ucov=0 121 ! vcov=0 122 ELSE 123 ! test non-tournant, onde se propageant en latitude 124 do j=1,jjp1 125 do i=1,iip1 126 ij=(j-1)*iip1+i 127 ps(ij) = 1e5*(1 + .1*exp(-100*(1+sin(rlatu(j)))**2) ) 128 enddo 129 enddo 130 131 ! rho = preff/(cpp*teta) 132 teta = .01*preff/cpp ! rho = 100 ; phi = ps/rho = 1e3 ; c=30 m/s = 2600 km/j = 23 degres / j 133 ucov=0. 134 vcov=0. 135 END IF 136 137 CALL pression ( ip1jmp1, ap, bp, ps, p ) 138 CALL massdair(p,masse) 139 140 END SUBROUTINE sw_case_williamson91_6 141 !-----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
