Changeset 2336 for LMDZ5/trunk
- Timestamp:
- Jul 31, 2015, 7:22:21 PM (9 years ago)
- Location:
- LMDZ5/trunk
- Files:
-
- 2 deleted
- 6 edited
- 16 moved
-
libf/dyn3d_common/caldyn0.F90 (modified) (2 diffs)
-
libf/dyn3d_common/convmas.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3d_common/convmas.F) (1 diff)
-
libf/dyn3d_common/enercin.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3d_common/enercin.F) (1 diff)
-
libf/dyn3d_common/flumass.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3d_common/flumass.F) (1 diff)
-
libf/dyn3d_common/massbar.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3d_common/massbar.F) (1 diff)
-
libf/dyn3d_common/massbarxy.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3d_common/massbarxy.F) (1 diff)
-
libf/dyn3d_common/tourpot.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3d_common/tourpot.F) (1 diff)
-
libf/dyn3d_common/vitvert.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3d_common/vitvert.F) (1 diff)
-
libf/dyn3dmem/ce0l.F90 (deleted)
-
libf/dyn3dmem/convmas1_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/convmas1_loc.F) (1 diff)
-
libf/dyn3dmem/convmas2_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/convmas2_loc.F) (1 diff)
-
libf/dyn3dmem/convmas_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/convmas_loc.F) (1 diff)
-
libf/dyn3dmem/enercin_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/enercin_loc.F) (1 diff)
-
libf/dyn3dmem/flumass_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/flumass_loc.F) (1 diff)
-
libf/dyn3dmem/massbar_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/massbar_loc.F) (1 diff)
-
libf/dyn3dmem/massbarxy_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/massbarxy_loc.F) (1 diff)
-
libf/dyn3dmem/tourpot_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/tourpot_loc.F) (1 diff)
-
libf/dyn3dmem/vitvert_loc.F90 (moved) (moved from LMDZ5/trunk/libf/dyn3dmem/vitvert_loc.F) (1 diff)
-
libf/dyn3dpar/ce0l.F90 (deleted)
-
libf/dynlonlat_phylonlat/phylmd/ce0l.F90 (modified) (3 diffs)
-
libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90 (modified) (23 diffs)
-
libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90 (modified) (22 diffs)
-
libf/phylmd/limit_netcdf.F90 (modified) (17 diffs)
-
makelmdz_fcm (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/dyn3d_common/caldyn0.F90
r2334 r2336 1 SUBROUTINE caldyn0(itau,ucov,vcov,teta,ps, pk,phis,phi,w,time)1 SUBROUTINE caldyn0(itau,ucov,vcov,teta,ps,masse,pk,phis,phi,w,pbaru,pbarv,time) 2 2 ! 3 3 !------------------------------------------------------------------------------- … … 15 15 !=============================================================================== 16 16 ! Arguments: 17 INTEGER, INTENT(IN) :: itau !--- 17 INTEGER, INTENT(IN) :: itau !--- TIME STEP INDEX 18 18 REAL, INTENT(IN) :: vcov (ip1jm ,llm) !--- V COVARIANT WIND 19 19 REAL, INTENT(IN) :: ucov (ip1jmp1 ,llm) !--- U COVARIANT WIND 20 20 REAL, INTENT(IN) :: teta (ip1jmp1 ,llm) !--- POTENTIAL TEMPERATURE 21 21 REAL, INTENT(IN) :: ps (ip1jmp1) !--- GROUND PRESSURE 22 REAL, INTENT(OUT) :: masse(ip1jmp1 ,llm) !--- MASS IN EACH CELL 22 23 REAL, INTENT(IN) :: pk (iip1,jjp1,llm) !--- PRESSURE 23 24 REAL, INTENT(IN) :: phis (ip1jmp1) !--- GROUND GEOPOTENTIAL 24 25 REAL, INTENT(IN) :: phi (ip1jmp1 ,llm) !--- 3D GEOPOTENTIAL 25 26 REAL, INTENT(OUT) :: w (ip1jmp1 ,llm) !--- VERTICAL WIND 27 REAL, INTENT(OUT) :: pbaru(ip1jmp1 ,llm) !--- U MASS FLUX 28 REAL, INTENT(OUT) :: pbarv(ip1jm ,llm) !--- V MASS FLUX 26 29 REAL, INTENT(IN) :: time !--- TIME 27 30 !=============================================================================== 28 31 ! Local variables: 29 REAL masse(ip1jmp1 ,llm) !--- MASS IN EACH CELL30 REAL pbaru(ip1jmp1 ,llm) !--- U MASS FLUX31 REAL pbarv(ip1jm ,llm) !--- V MASS FLUX32 32 REAL, DIMENSION(ip1jmp1,llmp1) :: p 33 33 REAL, DIMENSION(ip1jmp1,llm) :: ucont, massebx, ang, ecin, convm, bern -
LMDZ5/trunk/libf/dyn3d_common/convmas.F90
r2335 r2336 1 SUBROUTINE convmas (pbaru, pbarv, convm) 1 2 ! 2 ! $Header$ 3 ! 4 SUBROUTINE convmas (pbaru, pbarv, convm ) 5 c 6 IMPLICIT NONE 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute mass flux convergence at p levels. 7 IMPLICIT NONE 8 include "dimensions.h" 9 include "paramet.h" 10 include "comgeom.h" 11 include "logic.h" 12 !=============================================================================== 13 ! Arguments: 14 REAL, INTENT(IN) :: pbaru(ip1jmp1,llm) 15 REAL, INTENT(IN) :: pbarv(ip1jm ,llm) 16 REAL, INTENT(OUT) :: convm(ip1jmp1,llm) 17 !=============================================================================== 18 ! Method used: Computation from top to bottom. 19 ! Mass convergence at level llm is equal to zero and is not stored in convm. 20 !=============================================================================== 21 ! Local variables: 22 INTEGER :: l 23 !=============================================================================== 7 24 8 c======================================================================= 9 c 10 c Auteurs: P. Le Van , F. Hourdin . 11 c ------- 12 c 13 c Objet: 14 c ------ 15 c 16 c ******************************************************************** 17 c .... calcul de la convergence du flux de masse aux niveaux p ... 18 c ******************************************************************** 19 c 20 c 21 c pbaru et pbarv sont des arguments d'entree pour le s-pg .... 22 c ..... convm est un argument de sortie pour le s-pg .... 23 c 24 c le calcul se fait de haut en bas, 25 c la convergence de masse au niveau p(llm+1) est egale a 0. et 26 c n'est pas stockee dans le tableau convm . 27 c 28 c 29 c======================================================================= 30 c 31 c Declarations: 32 c ------------- 25 !--- Computation of - (d(pbaru)/dx + d(pbarv)/dy ) 26 CALL convflu( pbaru, pbarv, llm, convm ) 33 27 34 #include "dimensions.h" 35 #include "paramet.h" 36 #include "comvert.h" 37 #include "logic.h" 28 !--- Filter 29 CALL filtreg( convm, jjp1, llm, 2, 2, .TRUE., 1 ) 38 30 39 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ),convm( ip1jmp1,llm ) 40 INTEGER l,ij 31 !--- Mass convergence is integrated from top to bottom 32 DO l=llmm1,1,-1 33 convm(:,l) = convm(:,l) + convm(:,l+1) 34 END DO 41 35 42 43 c----------------------------------------------------------------------- 44 c .... calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ...... 45 46 CALL convflu( pbaru, pbarv, llm, convm ) 47 48 c----------------------------------------------------------------------- 49 c filtrage: 50 c --------- 51 52 CALL filtreg( convm, jjp1, llm, 2, 2, .true., 1 ) 53 54 c integration de la convergence de masse de haut en bas ...... 55 56 DO l = llmm1, 1, -1 57 DO ij = 1, ip1jmp1 58 convm(ij,l) = convm(ij,l) + convm(ij,l+1) 59 ENDDO 60 ENDDO 61 c 62 RETURN 63 END 36 END SUBROUTINE convmas -
LMDZ5/trunk/libf/dyn3d_common/enercin.F90
r2335 r2336 1 SUBROUTINE enercin ( vcov, ucov, vcont, ucont, ecin ) 1 2 ! 2 ! $Header$ 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute kinetic energy at sigma levels. 7 IMPLICIT NONE 8 include "dimensions.h" 9 include "paramet.h" 10 include "comgeom.h" 11 !=============================================================================== 12 ! Arguments: 13 REAL, INTENT(IN) :: vcov (ip1jm, llm) 14 REAL, INTENT(IN) :: ucov (ip1jmp1,llm) 15 REAL, INTENT(IN) :: vcont (ip1jm, llm) 16 REAL, INTENT(IN) :: ucont (ip1jmp1,llm) 17 REAL, INTENT(OUT) :: ecin (ip1jmp1,llm) 18 !=============================================================================== 19 ! Notes: 20 ! . V 21 ! i,j-1 3 22 ! 4 SUBROUTINE enercin ( vcov, ucov, vcont, ucont, ecin ) 5 IMPLICIT NONE 23 ! alpha4 . . alpha1 24 ! 25 ! 26 ! U . . P . U 27 ! i-1,j i,j i,j 28 ! 29 ! alpha3 . . alpha2 30 ! 31 ! 32 ! . V 33 ! i,j 34 ! 35 ! Kinetic energy at scalar point P(i,j) (excluding poles) is: 36 ! Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 ) + 37 ! 0.5 * U(i ,j)**2 *( alpha1 + alpha2 ) + 38 ! 0.5 * V(i,j-1)**2 *( alpha1 + alpha4 ) + 39 ! 0.5 * V(i, j)**2 *( alpha2 + alpha3 ) 40 !=============================================================================== 41 ! Local variables: 42 INTEGER :: l, ij, i 43 REAL :: ecinni(iip1), ecinsi(iip1), ecinpn, ecinps 44 !=============================================================================== 45 DO l=1,llm 46 DO ij = iip2, ip1jm -1 47 ecin(ij+1,l)=0.5*(ucov(ij ,l)*ucont(ij ,l)*alpha3p4(ij +1) & 48 + ucov(ij+1 ,l)*ucont(ij+1 ,l)*alpha1p2(ij +1) & 49 + vcov(ij-iim,l)*vcont(ij-iim,l)*alpha1p4(ij +1) & 50 + vcov(ij+1 ,l)*vcont(ij+1 ,l)*alpha2p3(ij +1) ) 51 END DO 52 !--- Correction: ecin(1,j,l)= ecin(iip1,j,l) 53 DO ij=iip2,ip1jm,iip1; ecin(ij,l) = ecin(ij+iim,l); END DO 6 54 7 c======================================================================= 8 c 9 c Auteur: P. Le Van 10 c ------- 11 c 12 c Objet: 13 c ------ 14 c 15 c ********************************************************************* 16 c .. calcul de l'energie cinetique aux niveaux s ...... 17 c ********************************************************************* 18 c vcov, vcont, ucov et ucont sont des arguments d'entree pour le s-pg . 19 c ecin est un argument de sortie pour le s-pg 20 c 21 c======================================================================= 55 !--- North pole 56 DO i=1,iim 57 ecinni(i) = vcov(i,l)*vcont(i,l)*aire(i) 58 END DO 59 ecinpn = 0.5*SUM(ecinni(1:iim))/apoln 60 DO ij=1,iip1; ecin(ij,l)=ecinpn; END DO 22 61 23 #include "dimensions.h" 24 #include "paramet.h" 25 #include "comgeom.h" 62 !--- South pole 63 DO i=1,iim 64 ecinsi(i) = vcov(i+ip1jmi1,l)*vcont(i+ip1jmi1,l)*aire(i+ip1jm) 65 END DO 66 ecinps = 0.5*SUM(ecinsi(1:iim))/apols 67 DO ij=1,iip1; ecin(ij+ip1jm,l)=ecinps; END DO 68 END DO 26 69 27 REAL vcov( ip1jm,llm ),vcont( ip1jm,llm ), 28 * ucov( ip1jmp1,llm ),ucont( ip1jmp1,llm ),ecin( ip1jmp1,llm ) 70 END SUBROUTINE enercin 29 71 30 REAL ecinni( iip1 ),ecinsi( iip1 )31 32 REAL ecinpn, ecinps33 INTEGER l,ij,i34 35 REAL SSUM36 37 38 39 c . V40 c i,j-141 42 c alpha4 . . alpha143 44 45 c U . . P . U46 c i-1,j i,j i,j47 48 c alpha3 . . alpha249 50 51 c . V52 c i,j53 54 c55 c L'energie cinetique au point scalaire P(i,j) ,autre que les poles, est :56 c Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 ) +57 c 0.5 * U(i ,j)**2 *( alpha1 + alpha2 ) +58 c 0.5 * V(i,j-1)**2 *( alpha1 + alpha4 ) +59 c 0.5 * V(i, j)**2 *( alpha2 + alpha3 )60 61 62 DO 5 l = 1,llm63 64 DO 1 ij = iip2, ip1jm -165 ecin( ij+1, l ) = 0.5 *66 * ( ucov( ij ,l ) * ucont( ij ,l ) * alpha3p4( ij +1 ) +67 * ucov( ij+1 ,l ) * ucont( ij+1 ,l ) * alpha1p2( ij +1 ) +68 * vcov(ij-iim,l ) * vcont(ij-iim,l ) * alpha1p4( ij +1 ) +69 * vcov( ij+ 1,l ) * vcont( ij+ 1,l ) * alpha2p3( ij +1 ) )70 1 CONTINUE71 72 c ... correction pour ecin(1,j,l) ....73 c ... ecin(1,j,l)= ecin(iip1,j,l) ...74 75 CDIR$ IVDEP76 DO 2 ij = iip2, ip1jm, iip177 ecin( ij,l ) = ecin( ij + iim, l )78 2 CONTINUE79 80 c calcul aux poles .......81 82 83 DO 3 i = 1, iim84 ecinni(i) = vcov( i , l) * vcont( i ,l) * aire( i )85 ecinsi(i) = vcov(i+ip1jmi1,l) * vcont(i+ip1jmi1,l) * aire(i+ip1jm)86 3 CONTINUE87 88 ecinpn = 0.5 * SSUM( iim,ecinni,1 ) / apoln89 ecinps = 0.5 * SSUM( iim,ecinsi,1 ) / apols90 91 DO 4 ij = 1,iip192 ecin( ij , l ) = ecinpn93 ecin( ij+ ip1jm, l ) = ecinps94 4 CONTINUE95 96 5 CONTINUE97 RETURN98 END -
LMDZ5/trunk/libf/dyn3d_common/flumass.F90
r2335 r2336 1 SUBROUTINE flumass (massebx,masseby, vcont, ucont, pbaru, pbarv ) 1 2 ! 2 ! $Header$ 3 ! 4 SUBROUTINE flumass (massebx,masseby, vcont, ucont, pbaru, pbarv ) 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute mass flux at s levels. 7 IMPLICIT NONE 8 include "dimensions.h" 9 include "paramet.h" 10 include "comgeom.h" 11 !=============================================================================== 12 ! Arguments: 13 REAL, INTENT(IN) :: massebx(ip1jmp1,llm) 14 REAL, INTENT(IN) :: masseby(ip1jm ,llm) 15 REAL, INTENT(IN) :: vcont (ip1jm ,llm) 16 REAL, INTENT(IN) :: ucont (ip1jmp1,llm) 17 REAL, INTENT(OUT) :: pbaru (ip1jmp1,llm) 18 REAL, INTENT(OUT) :: pbarv (ip1jm ,llm) 19 !=============================================================================== 20 ! Method used: A 2 equations system is solved. 21 ! * 1st one describes divergence computation at pole point nr. i (i=1 to im): 22 ! (0.5*(pbaru(i)-pbaru(i-1))-pbarv(i))/aire(i) = - SUM(pbarv(n))/aire pole 23 ! * 2nd one specifies that mean mass flux at pole is equal to 0: 24 ! SUM(pbaru(n)*local_area(n))=0 25 ! This way, we determine additive constant common to pbary elements representing 26 ! pbaru(0,j,l) in divergence computation equation for point i=1. (i=1 to im) 27 !=============================================================================== 28 ! Local variables: 29 REAL :: sairen, saireun, ctn, ctn0, apbarun(iip1) 30 REAL :: saires, saireus, cts, cts0, apbarus(iip1) 31 INTEGER :: l, i 32 !=============================================================================== 33 DO l=1,llm 34 pbaru(iip2:ip1jm,l)=massebx(iip2:ip1jm,l)*ucont(iip2:ip1jm,l) 35 pbarv( 1:ip1jm,l)=masseby( 1:ip1jm,l)*vcont( 1:ip1jm,l) 36 END DO 5 37 6 IMPLICIT NONE 38 !--- NORTH POLE 39 sairen =SUM(aire (1:iim)) 40 saireun=SUM(aireu(1:iim)) 41 DO l = 1,llm 42 ctn=SUM(pbarv(1:iim,l))/sairen 43 pbaru(1,l)= pbarv(1,l)-ctn*aire(1) 44 DO i=2,iim 45 pbaru(i,l)=pbaru(i-1,l)+pbarv(i,l)-ctn*aire(i) 46 END DO 47 DO i=1,iim 48 apbarun(i)=aireu(i)*pbaru(i,l) 49 END DO 50 ctn0 = -SUM(apbarun(1:iim))/saireun 51 DO i = 1,iim 52 pbaru(i,l)=2.*(pbaru(i,l)+ctn0) 53 END DO 54 pbaru(iip1,l)=pbaru(1,l) 55 END DO 7 56 8 c======================================================================= 9 c 10 c Auteurs: P. Le Van, F. Hourdin . 11 c ------- 12 c 13 c Objet: 14 c ------ 15 c 16 c ********************************************************************* 17 c .... calcul du flux de masse aux niveaux s ...... 18 c ********************************************************************* 19 c massebx,masseby,vcont et ucont sont des argum. d'entree pour le s-pg . 20 c pbaru et pbarv sont des argum.de sortie pour le s-pg . 21 c 22 c======================================================================= 57 !--- SOUTH POLE 58 saires =SUM(aire (ip1jm+1:ip1jmp1-1)) 59 saireus=SUM(aireu(ip1jm+1:ip1jmp1-1)) 60 DO l = 1,llm 61 cts=SUM(pbarv(ip1jmi1+1:ip1jm-1,l))/saires 62 pbaru(1+ip1jm,l)=-pbarv(1+ip1jmi1,l)+cts*aire(1+ip1jm) 63 DO i=2,iim 64 pbaru(i+ip1jm,l)=pbaru(i-1+ip1jm,l)-pbarv(i+ip1jmi1,l)+cts*aire(i+ip1jm) 65 END DO 66 DO i=1,iim 67 apbarus(i)=aireu(i+ip1jm)*pbaru(i+ip1jm,l) 68 END DO 69 cts0 = -SUM(apbarus(1:iim))/saireus 70 DO i = 1,iim 71 pbaru(i+ip1jm,l)=2.*(pbaru(i+ip1jm,l)+cts0) 72 END DO 73 pbaru(ip1jmp1,l)=pbaru(1+ip1jm,l) 74 END DO 23 75 24 25 #include "dimensions.h" 26 #include "paramet.h" 27 #include "comgeom.h" 28 29 REAL massebx( ip1jmp1,llm ),masseby( ip1jm,llm ) , 30 * vcont( ip1jm,llm ),ucont( ip1jmp1,llm ),pbaru( ip1jmp1,llm ), 31 * pbarv( ip1jm,llm ) 32 33 REAL apbarun( iip1 ),apbarus( iip1 ) 34 35 REAL sairen,saireun,saires,saireus,ctn,cts,ctn0,cts0 36 INTEGER l,ij,i 37 38 REAL SSUM 39 40 41 DO 5 l = 1,llm 42 43 DO 1 ij = iip2,ip1jm 44 pbaru( ij,l ) = massebx( ij,l ) * ucont( ij,l ) 45 1 CONTINUE 46 47 DO 3 ij = 1,ip1jm 48 pbarv( ij,l ) = masseby( ij,l ) * vcont( ij,l ) 49 3 CONTINUE 50 51 5 CONTINUE 52 53 c ................................................................ 54 c calcul de la composante du flux de masse en x aux poles ....... 55 c ................................................................ 56 c par la resolution d'1 systeme de 2 equations . 57 58 c la premiere equat.decrivant le calcul de la divergence en 1 point i 59 c du pole,ce calcul etant itere de i=1 a i=im . 60 c c.a.d , 61 c ( ( 0.5*pbaru(i)-0.5*pbaru(i-1) - pbarv(i))/aire(i) = 62 c - somme de ( pbarv(n) )/aire pole 63 64 c l'autre equat.specifiant que la moyenne du flux de masse au pole est =0. 65 c c.a.d somme de pbaru(n)*aire locale(n) = 0. 66 67 c on en revient ainsi a determiner la constante additive commune aux pbaru 68 c qui representait pbaru(0,j,l) dans l'equat.du calcul de la diverg.au pt 69 c i=1 . 70 c i variant de 1 a im 71 c n variant de 1 a im 72 73 sairen = SSUM( iim, aire( 1 ), 1 ) 74 saireun= SSUM( iim, aireu( 1 ), 1 ) 75 saires = SSUM( iim, aire( ip1jm+1 ), 1 ) 76 saireus= SSUM( iim, aireu( ip1jm+1 ), 1 ) 77 78 DO 20 l = 1,llm 79 80 ctn = SSUM( iim, pbarv( 1 ,l), 1 )/ sairen 81 cts = SSUM( iim, pbarv(ip1jmi1+ 1,l), 1 )/ saires 82 83 pbaru( 1 ,l )= pbarv( 1 ,l ) - ctn * aire( 1 ) 84 pbaru( ip1jm+1,l )= - pbarv( ip1jmi1+1,l ) + cts * aire( ip1jm+1 ) 85 86 DO 11 i = 2,iim 87 pbaru( i ,l ) = pbaru( i - 1 ,l ) + 88 * pbarv( i ,l ) - ctn * aire( i ) 89 90 pbaru( i+ ip1jm,l ) = pbaru( i+ ip1jm-1,l ) - 91 * pbarv( i+ ip1jmi1,l ) + cts * aire(i+ ip1jm) 92 11 CONTINUE 93 DO 12 i = 1,iim 94 apbarun(i) = aireu( i ) * pbaru( i , l) 95 apbarus(i) = aireu(i +ip1jm) * pbaru(i +ip1jm, l) 96 12 CONTINUE 97 ctn0 = -SSUM( iim,apbarun,1 )/saireun 98 cts0 = -SSUM( iim,apbarus,1 )/saireus 99 DO 14 i = 1,iim 100 pbaru( i , l) = 2. * ( pbaru( i , l) + ctn0 ) 101 pbaru(i+ ip1jm, l) = 2. * ( pbaru(i +ip1jm, l) + cts0 ) 102 14 CONTINUE 103 104 pbaru( iip1 ,l ) = pbaru( 1 ,l ) 105 pbaru( ip1jmp1,l ) = pbaru( ip1jm +1,l ) 106 20 CONTINUE 107 108 RETURN 109 END 76 END SUBROUTINE flumass -
LMDZ5/trunk/libf/dyn3d_common/massbar.F90
r2335 r2336 1 SUBROUTINE massbar(masse,massebx,masseby) 1 2 ! 2 ! $Header$ 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute air mass mean along X and Y in each cell. 7 ! See iniconst for more details. 8 IMPLICIT NONE 9 include "dimensions.h" 10 include "paramet.h" 11 include "comgeom.h" 12 !=============================================================================== 13 ! Arguments: 14 REAL, INTENT(IN) :: masse (ip1jmp1,llm) 15 REAL, INTENT(OUT) :: massebx(ip1jmp1,llm) 16 REAL, INTENT(OUT) :: masseby(ip1jm ,llm) 17 !------------------------------------------------------------------------------- 18 ! Method used. Each scalar point is associated to 4 area coefficients: 19 ! * alpha1(i,j) at point ( i+1/4,j-1/4 ) 20 ! * alpha2(i,j) at point ( i+1/4,j+1/4 ) 21 ! * alpha3(i,j) at point ( i-1/4,j+1/4 ) 22 ! * alpha4(i,j) at point ( i-1/4,j-1/4 ) 23 ! where alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j) 3 24 ! 4 SUBROUTINE massbar( masse, massebx, masseby ) 5 IMPLICIT NONE 6 c 7 c ********************************************************************** 8 c 9 c Calcule les moyennes en x et y de la masse d'air dans chaque maille. 10 c ********************************************************************** 11 c Auteurs : P. Le Van , Fr. Hourdin . 12 c .......... 13 c 14 c .. masse est un argum. d'entree pour le s-pg ... 15 c .. massebx,masseby sont des argum. de sortie pour le s-pg ... 16 c 17 c 18 c IMPLICIT NONE 19 c 20 #include "dimensions.h" 21 #include "paramet.h" 22 #include "comconst.h" 23 #include "comgeom.h" 24 c 25 REAL masse( ip1jmp1,llm ), massebx( ip1jmp1,llm ) , 26 * masseby( ip1jm,llm ) 27 INTEGER ij,l 28 c 29 c 30 c Methode pour calculer massebx et masseby . 31 c ---------------------------------------- 32 c 33 c A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires 34 c alpha1(i,j) calcule au point ( i+1/4,j-1/4 ) 35 c alpha2(i,j) calcule au point ( i+1/4,j+1/4 ) 36 c alpha3(i,j) calcule au point ( i-1/4,j+1/4 ) 37 c alpha4(i,j) calcule au point ( i-1/4,j-1/4 ) 38 c 39 c Avec alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j) 40 c 41 c N.B . Pour plus de details, voir s-pg ... iniconst ... 42 c 43 c 44 c 45 c alpha4 . . alpha1 . alpha4 46 c (i,j) (i,j) (i+1,j) 47 c 48 c P . U . . P 49 c (i,j) (i,j) (i+1,j) 50 c 51 c alpha3 . . alpha2 .alpha3 52 c (i,j) (i,j) (i+1,j) 53 c 54 c V . Z . . V 55 c (i,j) 56 c 57 c alpha4 . . alpha1 .alpha4 58 c (i,j+1) (i,j+1) (i+1,j+1) 59 c 60 c P . U . . P 61 c (i,j+1) (i+1,j+1) 62 c 63 c 64 c 65 c On a : 66 c 67 c massebx(i,j) = masse(i ,j) * ( alpha1(i ,j) + alpha2(i,j)) + 68 c masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) ) 69 c localise au point ... U (i,j) ... 70 c 71 c masseby(i,j) = masse(i,j ) * ( alpha2(i,j ) + alpha3(i,j ) + 72 c masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1) 73 c localise au point ... V (i,j) ... 74 c 75 c 76 c======================================================================= 25 ! alpha4 . . alpha1 . alpha4 26 ! (i,j) (i,j) (i+1,j) 27 ! 28 ! P . U . . P 29 ! (i,j) (i,j) (i+1,j) 30 ! 31 ! alpha3 . . alpha2 .alpha3 32 ! (i,j) (i,j) (i+1,j) 33 ! 34 ! V . Z . . V 35 ! (i,j) 36 ! 37 ! alpha4 . . alpha1 .alpha4 38 ! (i,j+1) (i,j+1) (i+1,j+1) 39 ! 40 ! P . U . . P 41 ! (i,j+1) (i+1,j+1) 42 ! 43 ! 44 ! massebx(i,j) = masse(i ,j) * ( alpha1(i ,j) + alpha2(i,j)) + 45 ! masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) ) 46 ! localized at point ... U (i,j) ... 47 ! 48 ! masseby(i,j) = masse(i,j ) * ( alpha2(i,j ) + alpha3(i,j ) + 49 ! masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1) 50 ! localized at point ... V (i,j) ... 51 !=============================================================================== 52 ! Local variables: 53 INTEGER :: ij, l 54 !=============================================================================== 55 DO l=1,llm 56 DO ij=1,ip1jmp1-1 57 massebx(ij,l)=masse(ij,l)*alpha1p2(ij)+masse(ij+1 ,l)*alpha3p4(ij+1) 58 END DO 59 DO ij=iip1,ip1jmp1,iip1; massebx(ij,l)=massebx(ij-iim,l); END DO 60 DO ij=1,ip1jm 61 masseby(ij,l)=masse(ij,l)*alpha2p3(ij)+masse(ij+iip1,l)*alpha1p4(ij+iip1) 62 END DO 63 END DO 77 64 78 DO 100 l = 1 , llm 79 c 80 DO ij = 1, ip1jmp1 - 1 81 massebx(ij,l) = masse( ij, l) * alpha1p2( ij ) + 82 * masse(ij+1, l) * alpha3p4(ij+1 ) 83 ENDDO 65 END SUBROUTINE massbar 84 66 85 c .... correction pour massebx( iip1,j) .....86 c ... massebx(iip1,j)= massebx(1,j) ...87 c88 CDIR$ IVDEP89 DO ij = iip1, ip1jmp1, iip190 massebx( ij,l ) = massebx( ij - iim,l )91 ENDDO92 93 94 DO ij = 1,ip1jm95 masseby( ij,l ) = masse( ij , l ) * alpha2p3( ij ) +96 * masse(ij+iip1, l ) * alpha1p4( ij+iip1 )97 ENDDO98 99 100 CONTINUE100 c101 RETURN102 END -
LMDZ5/trunk/libf/dyn3d_common/massbarxy.F90
r2335 r2336 1 SUBROUTINE massbarxy(masse,massebxy) 1 2 ! 2 ! $Header$ 3 ! 4 SUBROUTINE massbarxy( masse, massebxy ) 5 IMPLICIT NONE 6 c 7 c ********************************************************************** 8 c 9 c Calcule les moyennes en x et y de la masse d'air dans chaque maille. 10 c ********************************************************************** 11 c Auteurs : P. Le Van , Fr. Hourdin . 12 c .......... 13 c 14 c .. masse est un argum. d'entree pour le s-pg ... 15 c .. massebxy est un argum. de sortie pour le s-pg ... 16 c 17 c 18 c IMPLICIT NONE 19 c 20 #include "dimensions.h" 21 #include "paramet.h" 22 #include "comconst.h" 23 #include "comgeom.h" 24 c 25 REAL masse( ip1jmp1,llm ), massebxy( ip1jm,llm ) 26 INTEGER ij,l 27 c 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute air mass mean along X and Y in each cell. 7 ! See iniconst for more details. 8 IMPLICIT NONE 9 include "dimensions.h" 10 include "paramet.h" 11 include "comconst.h" 12 include "comgeom.h" 13 !=============================================================================== 14 ! Arguments: 15 REAL, INTENT(IN) :: masse (ip1jmp1,llm) 16 REAL, INTENT(OUT) :: massebxy(ip1jm ,llm) 17 !=============================================================================== 18 ! Local variables: 19 INTEGER :: ij, l 20 !=============================================================================== 21 DO l=1,llm 22 DO ij=1,ip1jm-1 23 massebxy(ij,l)=masse(ij ,l)*alpha2(ij ) + & 24 + masse(ij+1 ,l)*alpha3(ij+1 ) + & 25 + masse(ij+iip1,l)*alpha1(ij+iip1) + & 26 + masse(ij+iip2,l)*alpha4(ij+iip2) 27 END DO 28 DO ij=iip1,ip1jm,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO 29 END DO 28 30 29 DO 100 l = 1 , llm 30 c 31 DO 5 ij = 1, ip1jm - 1 32 massebxy( ij,l ) = masse( ij ,l ) * alpha2( ij ) + 33 + masse( ij+1 ,l ) * alpha3( ij+1 ) + 34 + masse( ij+iip1,l ) * alpha1( ij+iip1 ) + 35 + masse( ij+iip2,l ) * alpha4( ij+iip2 ) 36 5 CONTINUE 37 38 c .... correction pour massebxy( iip1,j ) ........ 39 40 CDIR$ IVDEP 41 42 DO 7 ij = iip1, ip1jm, iip1 43 massebxy( ij,l ) = massebxy( ij - iim,l ) 44 7 CONTINUE 45 46 100 CONTINUE 47 c 48 RETURN 49 END 31 END SUBROUTINE massbarxy -
LMDZ5/trunk/libf/dyn3d_common/tourpot.F90
r2335 r2336 1 SUBROUTINE tourpot ( vcov, ucov, massebxy, vorpot ) 1 2 ! 2 ! $Header$ 3 ! 4 SUBROUTINE tourpot ( vcov, ucov, massebxy, vorpot ) 5 IMPLICIT NONE 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute potential vorticity. 7 IMPLICIT NONE 8 include "dimensions.h" 9 include "paramet.h" 10 include "comgeom.h" 11 include "logic.h" 12 !=============================================================================== 13 ! Arguments: 14 REAL, INTENT(IN) :: vcov (ip1jm, llm) 15 REAL, INTENT(IN) :: ucov (ip1jmp1,llm) 16 REAL, INTENT(IN) :: massebxy(ip1jm, llm) 17 REAL, INTENT(OUT) :: vorpot (ip1jm, llm) 18 !=============================================================================== 19 ! Method used: 20 ! vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy 21 !=============================================================================== 22 ! Local variables: 23 INTEGER :: l, ij 24 REAL :: rot(ip1jm,llm) 25 !=============================================================================== 6 26 7 c======================================================================= 8 c 9 c Auteur: P. Le Van 10 c ------- 11 c 12 c Objet: 13 c ------ 14 c 15 c ******************************************************************* 16 c ......... calcul du tourbillon potentiel ......... 17 c ******************************************************************* 18 c 19 c vcov,ucov,fext et pbarxyfl sont des argum. d'entree pour le s-pg . 20 c vorpot est un argum.de sortie pour le s-pg . 21 c 22 c======================================================================= 27 !--- Wind vorticity ; correction: rot(iip1,j,l) = rot(1,j,l) 28 DO l=1,llm 29 DO ij=1,ip1jm-1 30 rot(ij,l)=vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l) 31 END DO 32 DO ij=iip1,ip1jm,iip1; rot(ij,l)=rot(ij-iim,l); END DO 33 END DO 23 34 24 #include "dimensions.h" 25 #include "paramet.h" 26 #include "comgeom.h" 27 #include "logic.h" 35 !--- Filter 36 CALL filtreg(rot,jjm,llm,2,1,.FALSE.,1) 28 37 29 REAL rot( ip1jm,llm ) 30 REAL vcov( ip1jm,llm ),ucov( ip1jmp1,llm ) 31 REAL massebxy( ip1jm,llm ),vorpot( ip1jm,llm ) 38 !--- Potential vorticity ; correction: rot(iip1,j,l) = rot(1,j,l) 39 DO l=1,llm 40 DO ij=1,ip1jm-1 41 vorpot(ij,l)=(rot(ij,l)+fext(ij))/massebxy(ij,l) 42 END DO 43 DO ij=iip1,ip1jm,iip1; vorpot(ij,l)=vorpot(ij-iim,l); END DO 44 END DO 32 45 33 INTEGER l, ij 34 35 36 37 38 c ... vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy .. 39 40 41 42 c ........ Calcul du rotationnel du vent V puis filtrage ........ 43 44 DO 5 l = 1,llm 45 46 DO 2 ij = 1, ip1jm - 1 47 rot( ij,l ) = vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l) 48 2 CONTINUE 49 50 c .... correction pour rot( iip1,j,l ) ..... 51 c .... rot(iip1,j,l) = rot(1,j,l) ..... 52 53 CDIR$ IVDEP 54 55 DO 3 ij = iip1, ip1jm, iip1 56 rot( ij,l ) = rot( ij -iim, l ) 57 3 CONTINUE 58 59 5 CONTINUE 60 61 62 CALL filtreg( rot, jjm, llm, 2, 1, .FALSE., 1 ) 63 64 65 DO 10 l = 1, llm 66 67 DO 6 ij = 1, ip1jm - 1 68 vorpot( ij,l ) = ( rot(ij,l) + fext(ij) ) / massebxy(ij,l) 69 6 CONTINUE 70 71 c ..... correction pour vorpot( iip1,j,l) ..... 72 c .... vorpot(iip1,j,l)= vorpot(1,j,l) .... 73 CDIR$ IVDEP 74 DO 8 ij = iip1, ip1jm, iip1 75 vorpot( ij,l ) = vorpot( ij -iim,l ) 76 8 CONTINUE 77 78 10 CONTINUE 79 80 RETURN 81 END 46 END SUBROUTINE tourpot -
LMDZ5/trunk/libf/dyn3d_common/vitvert.F90
r2335 r2336 1 SUBROUTINE vitvert (convm, w) 1 2 ! 2 ! $Header$ 3 ! 4 SUBROUTINE vitvert ( convm , w ) 5 c 6 IMPLICIT NONE 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute vertical speed at sigma levels. 7 IMPLICIT NONE 8 include "dimensions.h" 9 include "paramet.h" 10 include "comvert.h" 11 !=============================================================================== 12 ! Arguments: 13 REAL, INTENT(IN) :: convm(ip1jmp1,llm) 14 REAL, INTENT(OUT) :: w (ip1jmp1,llm) 15 !=============================================================================== 16 ! Notes: Vertical speed is oriented from bottom to top. 17 ! * At ground - level sigma(1): w(i,j,1) = 0. 18 ! * At top - level sigma(llm+1): w(i,j,l) = 0. (not stored in w) 19 !=============================================================================== 20 ! Local variables: 21 INTEGER :: l 22 !=============================================================================== 23 DO l=1,llmm1; w(:,l+1)=convm(:,l+1)-bp(l+1)*convm(:,1); END DO 24 w(:,1)=0. 7 25 8 c======================================================================= 9 c 10 c Auteurs: P. Le Van , F. Hourdin . 11 c ------- 12 c 13 c Objet: 14 c ------ 15 c 16 c ******************************************************************* 17 c .... calcul de la vitesse verticale aux niveaux sigma .... 18 c ******************************************************************* 19 c convm est un argument d'entree pour le s-pg ...... 20 c w est un argument de sortie pour le s-pg ...... 21 c 22 c la vitesse verticale est orientee de haut en bas . 23 c au sol, au niveau sigma(1), w(i,j,1) = 0. 24 c au sommet, au niveau sigma(llm+1) , la vit.verticale est aussi 25 c egale a 0. et n'est pas stockee dans le tableau w . 26 c 27 c 28 c======================================================================= 26 END SUBROUTINE vitvert 29 27 30 #include "dimensions.h"31 #include "paramet.h"32 #include "comvert.h"33 34 REAL w(ip1jmp1,llm),convm(ip1jmp1,llm)35 INTEGER l, ij36 37 38 39 DO 2 l = 1,llmm140 41 DO 1 ij = 1,ip1jmp142 w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )43 1 CONTINUE44 45 2 CONTINUE46 47 DO 5 ij = 1,ip1jmp148 w(ij,1) = 0.49 5 CONTINUE50 51 RETURN52 END -
LMDZ5/trunk/libf/dyn3dmem/convmas1_loc.F90
r2335 r2336 1 SUBROUTINE convmas1_loc (pbaru, pbarv, convm ) 2 c 3 USE parallel_lmdz 4 USE mod_filtreg_p 5 IMPLICIT NONE 1 SUBROUTINE convmas1_loc (pbaru, pbarv, convm) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute mass flux convergence at p levels. 7 ! Equivalent to convmas_loc if convmas2_loc is called after. 8 USE parallel_lmdz 9 USE mod_filtreg_p 10 IMPLICIT NONE 11 include "dimensions.h" 12 include "paramet.h" 13 include "comgeom.h" 14 include "logic.h" 15 !=============================================================================== 16 ! Arguments: 17 REAL, INTENT(IN) :: pbaru(ijb_u:ije_u,llm) 18 REAL, INTENT(IN) :: pbarv(ijb_v:ije_v,llm) 19 REAL, TARGET, INTENT(OUT) :: convm(ijb_u:ije_u,llm) 20 !=============================================================================== 21 ! Method used: Computation from top to bottom. 22 ! Mass convergence at level llm is equal to zero and is not stored in convm. 23 !=============================================================================== 24 ! Local variables: 25 INTEGER :: l, jjb, jje 26 !=============================================================================== 6 27 7 c======================================================================= 8 c 9 c Auteurs: P. Le Van , F. Hourdin . 10 c ------- 11 c 12 c Objet: 13 c ------ 14 c 15 c ******************************************************************** 16 c .... calcul de la convergence du flux de masse aux niveaux p ... 17 c ******************************************************************** 18 c 19 c 20 c pbaru et pbarv sont des arguments d'entree pour le s-pg .... 21 c ..... convm est un argument de sortie pour le s-pg .... 22 c 23 c le calcul se fait de haut en bas, 24 c la convergence de masse au niveau p(llm+1) est egale a 0. et 25 c n'est pas stockee dans le tableau convm . 26 c 27 c 28 c======================================================================= 29 c 30 c Declarations: 31 c ------------- 28 !--- Computation of - (d(pbaru)/dx + d(pbarv)/dy ) 29 CALL convflu_loc( pbaru, pbarv, llm, convm ) 32 30 33 #include "dimensions.h" 34 #include "paramet.h" 35 #include "comvert.h" 36 #include "logic.h" 31 !--- Filter 32 jjb=jj_begin 33 jje=jj_end+1 34 IF(pole_sud) jje=jj_end 35 CALL filtreg_p(convm,jjb_u,jje_u,jjb,jje,jjp1,llm,2,2,.TRUE.,1) 37 36 38 REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm ) 39 REAL, target :: convm( ijb_u:ije_u,llm ) 40 INTEGER l,ij 37 END SUBROUTINE convmas1_loc 41 38 42 INTEGER ijb,ije,jjb,jje43 44 45 c-----------------------------------------------------------------------46 c .... calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ......47 48 CALL convflu_loc( pbaru, pbarv, llm, convm )49 50 c-----------------------------------------------------------------------51 c filtrage:52 c ---------53 54 jjb=jj_begin55 jje=jj_end+156 if (pole_sud) jje=jj_end57 58 CALL filtreg_p( convm, jjb_u,jje_u,jjb, jje, jjp1, llm,59 & 2, 2, .true., 1 )60 61 c integration de la convergence de masse de haut en bas ......62 c63 RETURN64 END -
LMDZ5/trunk/libf/dyn3dmem/convmas2_loc.F90
r2335 r2336 1 SUBROUTINE convmas2_loc ( convm ) 2 c 3 USE parallel_lmdz 4 IMPLICIT NONE 1 SUBROUTINE convmas2_loc (convm) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute mass flux convergence at p levels. 7 ! Equivalent to convmas_loc if convmas1_loc is called before. 8 USE parallel_lmdz 9 IMPLICIT NONE 10 include "dimensions.h" 11 include "paramet.h" 12 include "comgeom.h" 13 include "logic.h" 14 !=============================================================================== 15 ! Arguments: 16 REAL, INTENT(INOUT) :: convm(ijb_u:ije_u,llm) 17 !=============================================================================== 18 ! Method used: Computation from top to bottom. 19 ! Mass convergence at level llm is equal to zero and is not stored in convm. 20 !=============================================================================== 21 ! Local variables: 22 INTEGER :: l, ijb, ije 23 !=============================================================================== 5 24 6 c======================================================================= 7 c 8 c Auteurs: P. Le Van , F. Hourdin . 9 c ------- 10 c 11 c Objet: 12 c ------ 13 c 14 c ******************************************************************** 15 c .... calcul de la convergence du flux de masse aux niveaux p ... 16 c ******************************************************************** 17 c 18 c 19 c pbaru et pbarv sont des arguments d'entree pour le s-pg .... 20 c ..... convm est un argument de sortie pour le s-pg .... 21 c 22 c le calcul se fait de haut en bas, 23 c la convergence de masse au niveau p(llm+1) est egale a 0. et 24 c n'est pas stockee dans le tableau convm . 25 c 26 c 27 c======================================================================= 28 c 29 c Declarations: 30 c ------------- 25 !$OMP MASTER 26 !--- Mass convergence is integrated from top to bottom 27 ijb=ij_begin 28 ije=ij_end+iip1 29 IF(pole_sud) ije=ij_end 30 DO l=llmm1,1,-1 31 convm(ijb:ije,l) = convm(ijb:ije,l) + convm(ijb:ije,l+1) 32 END DO 33 !$OMP END MASTER 31 34 32 #include "dimensions.h" 33 #include "paramet.h" 34 #include "comvert.h" 35 #include "logic.h" 35 END SUBROUTINE convmas2_loc 36 36 37 REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm )38 REAL :: convm( ijb_u:ije_u,llm )39 INTEGER l,ij40 INTEGER ijb,ije,jjb,jje41 42 c$OMP MASTER43 c integration de la convergence de masse de haut en bas ......44 ijb=ij_begin45 ije=ij_end+iip146 if (pole_sud) ije=ij_end47 48 DO l = llmm1, 1, -149 DO ij = ijb, ije50 convm(ij,l) = convm(ij,l) + convm(ij,l+1)51 ENDDO52 ENDDO53 c54 c$OMP END MASTER55 RETURN56 END -
LMDZ5/trunk/libf/dyn3dmem/convmas_loc.F90
r2335 r2336 1 SUBROUTINE convmas_loc (pbaru, pbarv, convm ) 2 c 3 USE parallel_lmdz 4 USE mod_filtreg_p 5 IMPLICIT NONE 1 SUBROUTINE convmas_loc (pbaru, pbarv, convm) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute mass flux convergence at p levels. 7 USE parallel_lmdz 8 USE mod_filtreg_p 9 IMPLICIT NONE 10 include "dimensions.h" 11 include "paramet.h" 12 include "comgeom.h" 13 include "logic.h" 14 !=============================================================================== 15 ! Arguments: 16 REAL, INTENT(IN) :: pbaru(ijb_u:ije_u,llm) 17 REAL, INTENT(IN) :: pbarv(ijb_v:ije_v,llm) 18 REAL, INTENT(OUT) :: convm(ijb_u:ije_u,llm) 19 !=============================================================================== 20 ! Method used: Computation from top to bottom. 21 ! Mass convergence at level llm is equal to zero and is not stored in convm. 22 !=============================================================================== 23 ! Local variables: 24 INTEGER :: l, ijb, ije, jjb, jje 25 !=============================================================================== 6 26 7 c======================================================================= 8 c 9 c Auteurs: P. Le Van , F. Hourdin . 10 c ------- 11 c 12 c Objet: 13 c ------ 14 c 15 c ******************************************************************** 16 c .... calcul de la convergence du flux de masse aux niveaux p ... 17 c ******************************************************************** 18 c 19 c 20 c pbaru et pbarv sont des arguments d'entree pour le s-pg .... 21 c ..... convm est un argument de sortie pour le s-pg .... 22 c 23 c le calcul se fait de haut en bas, 24 c la convergence de masse au niveau p(llm+1) est egale a 0. et 25 c n'est pas stockee dans le tableau convm . 26 c 27 c 28 c======================================================================= 29 c 30 c Declarations: 31 c ------------- 27 !--- Computation of - (d(pbaru)/dx + d(pbarv)/dy ) 28 CALL convflu_loc( pbaru, pbarv, llm, convm ) 32 29 33 #include "dimensions.h" 34 #include "paramet.h" 35 #include "comvert.h" 36 #include "logic.h" 30 !--- Filter 31 jjb=jj_begin 32 jje=jj_end+1 33 IF(pole_sud) jje=jj_end 34 CALL filtreg_p(convm,jjb_u,jje_u,jjb,jje,jjp1,llm,2,2,.TRUE.,1) 37 35 38 REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm ) 39 REAL, target :: convm( ijb_u:ije_u,llm ) 40 INTEGER l,ij 41 42 INTEGER ijb,ije,jjb,jje 43 44 45 c----------------------------------------------------------------------- 46 c .... calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ...... 47 48 CALL convflu_loc( pbaru, pbarv, llm, convm ) 49 50 c----------------------------------------------------------------------- 51 c filtrage: 52 c --------- 53 54 jjb=jj_begin 55 jje=jj_end+1 56 if (pole_sud) jje=jj_end 57 58 CALL filtreg_p(convm, jjb_u, jje_u,jjb, jje, jjp1, llm, 59 & 2, 2, .true., 1 ) 60 61 c integration de la convergence de masse de haut en bas ...... 36 !--- Mass convergence is integrated from top to bottom 62 37 !$OMP BARRIER 63 38 !$OMP MASTER 64 ijb=ij_begin 65 ije=ij_end+iip1 66 if (pole_sud) ije=ij_end 67 68 DO l = llmm1, 1, -1 69 DO ij = ijb, ije 70 convm(ij,l) = convm(ij,l) + convm(ij,l+1) 71 ENDDO 72 ENDDO 73 c 39 ijb=ij_begin 40 ije=ij_end+iip1 41 IF(pole_sud) ije=ij_end 42 DO l=llmm1,1,-1 43 convm(ijb:ije,l) = convm(ijb:ije,l) + convm(ijb:ije,l+1) 44 END DO 74 45 !$OMP END MASTER 75 46 !$OMP BARRIER 76 RETURN 77 END 47 48 END SUBROUTINE convmas_loc 49 -
LMDZ5/trunk/libf/dyn3dmem/enercin_loc.F90
r2335 r2336 1 SUBROUTINE enercin_loc ( vcov, ucov, vcont, ucont, ecin ) 2 USE parallel_lmdz 3 IMPLICIT NONE 1 SUBROUTINE enercin_loc ( vcov, ucov, vcont, ucont, ecin ) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute kinetic energy at sigma levels. 7 USE parallel_lmdz 8 IMPLICIT NONE 9 include "dimensions.h" 10 include "paramet.h" 11 include "comgeom.h" 12 !=============================================================================== 13 ! Arguments: 14 REAL, INTENT(IN) :: vcov (ijb_v:ije_v,llm) 15 REAL, INTENT(IN) :: ucov (ijb_u:ije_u,llm) 16 REAL, INTENT(IN) :: vcont (ijb_v:ije_v,llm) 17 REAL, INTENT(IN) :: ucont (ijb_u:ije_u,llm) 18 REAL, INTENT(OUT) :: ecin (ijb_u:ije_u,llm) 19 !=============================================================================== 20 ! Notes: 21 ! . V 22 ! i,j-1 23 ! 24 ! alpha4 . . alpha1 25 ! 26 ! 27 ! U . . P . U 28 ! i-1,j i,j i,j 29 ! 30 ! alpha3 . . alpha2 31 ! 32 ! 33 ! . V 34 ! i,j 35 ! 36 ! Kinetic energy at scalar point P(i,j) (excluding poles) is: 37 ! Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 ) + 38 ! 0.5 * U(i ,j)**2 *( alpha1 + alpha2 ) + 39 ! 0.5 * V(i,j-1)**2 *( alpha1 + alpha4 ) + 40 ! 0.5 * V(i, j)**2 *( alpha2 + alpha3 ) 41 !=============================================================================== 42 ! Local variables: 43 INTEGER :: l, ij, i, ijb, ije 44 REAL :: ecinni(iim), ecinsi(iim), ecinpn, ecinps 45 !=============================================================================== 46 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 47 DO l=1,llm 4 48 5 c======================================================================= 6 c 7 c Auteur: P. Le Van 8 c ------- 9 c 10 c Objet: 11 c ------ 12 c 13 c ********************************************************************* 14 c .. calcul de l'energie cinetique aux niveaux s ...... 15 c ********************************************************************* 16 c vcov, vcont, ucov et ucont sont des arguments d'entree pour le s-pg . 17 c ecin est un argument de sortie pour le s-pg 18 c 19 c======================================================================= 49 ijb=ij_begin 50 ije=ij_end+iip1 20 51 21 #include "dimensions.h" 22 #include "paramet.h" 23 #include "comgeom.h" 52 IF(pole_nord) ijb=ij_begin+iip1 53 IF(pole_sud) ije=ij_end-iip1 24 54 25 REAL vcov( ijb_v:ije_v,llm ),vcont( ijb_v:ije_v,llm ) 26 REAL ucov( ijb_u:ije_u,llm ),ucont( ijb_u:ije_u,llm ) 27 REAL ecin( ijb_u:ije_u,llm ) 55 DO ij = ijb,ije-1 56 ecin(ij+1,l)=0.5*(ucov(ij ,l)*ucont(ij ,l)*alpha3p4(ij +1) & 57 + ucov(ij+1 ,l)*ucont(ij+1 ,l)*alpha1p2(ij +1) & 58 + vcov(ij-iim,l)*vcont(ij-iim,l)*alpha1p4(ij +1) & 59 + vcov(ij+1 ,l)*vcont(ij+1 ,l)*alpha2p3(ij +1) ) 60 END DO 28 61 29 REAL ecinni( iip1 ),ecinsi( iip1 ) 62 !--- Correction: ecin(1,j,l)= ecin(iip1,j,l) 63 DO ij=ijb,ije,iip1; ecin(ij,l) = ecin(ij+iim,l); END DO 30 64 31 REAL ecinpn, ecinps 32 INTEGER l,ij,i,ijb,ije 65 !--- North pole 66 IF(pole_nord) THEN 67 ecinni(:) = vcov(1:iim,l)*vcont(1:iim,l)*aire(1:iim) 68 ecinpn = 0.5*SUM(ecinni)/apoln 69 ecin(1:iip1,l)=ecinpn 70 END IF 33 71 34 EXTERNAL SSUM 35 REAL SSUM 72 !--- South pole 73 IF(pole_sud) THEN 74 DO i=1,iim 75 ecinsi(i) = vcov(i+ip1jmi1,l)*vcont(i+ip1jmi1,l)*aire(i+ip1jm) 76 END DO 77 ecinps = 0.5*SUM(ecinsi)/apols 78 ecin(1+ip1jm:ip1jmp1,l)=ecinps 79 END IF 80 END DO 81 !$OMP END DO NOWAIT 36 82 83 END SUBROUTINE enercin_loc 37 84 38 39 c . V40 c i,j-141 42 c alpha4 . . alpha143 44 45 c U . . P . U46 c i-1,j i,j i,j47 48 c alpha3 . . alpha249 50 51 c . V52 c i,j53 54 c55 c L'energie cinetique au point scalaire P(i,j) ,autre que les poles, est :56 c Ecin = 0.5 * U(i-1,j)**2 *( alpha3 + alpha4 ) +57 c 0.5 * U(i ,j)**2 *( alpha1 + alpha2 ) +58 c 0.5 * V(i,j-1)**2 *( alpha1 + alpha4 ) +59 c 0.5 * V(i, j)**2 *( alpha2 + alpha3 )60 61 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)62 DO 5 l = 1,llm63 64 ijb=ij_begin65 ije=ij_end+iip166 67 IF (pole_nord) ijb=ij_begin+iip168 IF (pole_sud) ije=ij_end-iip169 70 DO 1 ij = ijb, ije -171 ecin( ij+1, l ) = 0.5 *72 * ( ucov( ij ,l ) * ucont( ij ,l ) * alpha3p4( ij +1 ) +73 * ucov( ij+1 ,l ) * ucont( ij+1 ,l ) * alpha1p2( ij +1 ) +74 * vcov(ij-iim,l ) * vcont(ij-iim,l ) * alpha1p4( ij +1 ) +75 * vcov( ij+ 1,l ) * vcont( ij+ 1,l ) * alpha2p3( ij +1 ) )76 1 CONTINUE77 78 c ... correction pour ecin(1,j,l) ....79 c ... ecin(1,j,l)= ecin(iip1,j,l) ...80 81 CDIR$ IVDEP82 DO 2 ij = ijb, ije, iip183 ecin( ij,l ) = ecin( ij + iim, l )84 2 CONTINUE85 86 c calcul aux poles .......87 88 IF (pole_nord) THEN89 90 DO i = 1, iim91 ecinni(i) = vcov( i , l) *92 * vcont( i ,l) * aire( i )93 ENDDO94 95 ecinpn = 0.5 * SSUM( iim,ecinni,1 ) / apoln96 97 DO ij = 1,iip198 ecin( ij , l ) = ecinpn99 ENDDO100 101 ENDIF102 103 IF (pole_sud) THEN104 105 DO i = 1, iim106 ecinsi(i) = vcov(i+ip1jmi1,l)*107 * vcont(i+ip1jmi1,l) * aire(i+ip1jm)108 ENDDO109 110 ecinps = 0.5 * SSUM( iim,ecinsi,1 ) / apols111 112 DO ij = 1,iip1113 ecin( ij+ ip1jm, l ) = ecinps114 ENDDO115 116 ENDIF117 118 119 5 CONTINUE120 c$OMP END DO NOWAIT121 RETURN122 END -
LMDZ5/trunk/libf/dyn3dmem/flumass_loc.F90
r2335 r2336 1 SUBROUTINE flumass_loc(massebx,masseby,vcont,ucont,pbaru,pbarv) 2 USE parallel_lmdz 3 IMPLICIT NONE 1 SUBROUTINE flumass_loc(massebx,masseby, vcont, ucont, pbaru, pbarv ) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute mass flux at s levels. 7 USE parallel_lmdz 8 IMPLICIT NONE 9 include "dimensions.h" 10 include "paramet.h" 11 include "comgeom.h" 12 !=============================================================================== 13 ! Arguments: 14 REAL, INTENT(IN) :: massebx(ijb_u:ije_u,llm) 15 REAL, INTENT(IN) :: masseby(ijb_v:ije_v,llm) 16 REAL, INTENT(IN) :: vcont (ijb_v:ije_v,llm) 17 REAL, INTENT(IN) :: ucont (ijb_u:ije_u,llm) 18 REAL, INTENT(OUT) :: pbaru (ijb_u:ije_u,llm) 19 REAL, INTENT(OUT) :: pbarv (ijb_v:ije_v,llm) 20 !=============================================================================== 21 ! Method used: A 2 equations system is solved. 22 ! * 1st one describes divergence computation at pole point nr. i (i=1 to im): 23 ! (0.5*(pbaru(i)-pbaru(i-1))-pbarv(i))/aire(i) = - SUM(pbarv(n))/aire pole 24 ! * 2nd one specifies that mean mass flux at pole is equal to 0: 25 ! SUM(pbaru(n)*local_area(n))=0 26 ! This way, we determine additive constant common to pbary elements representing 27 ! pbaru(0,j,l) in divergence computation equation for point i=1. (i=1 to im) 28 !=============================================================================== 29 ! Local variables: 30 REAL :: sairen, saireun, ctn, ctn0, apbarun(iim) 31 REAL :: saires, saireus, cts, cts0, apbarus(iim) 32 INTEGER :: l, i, ij, ijb, ije 33 !=============================================================================== 34 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 35 DO l=1,llm 4 36 5 c======================================================================= 6 c 7 c Auteurs: P. Le Van, F. Hourdin . 8 c ------- 9 c 10 c Objet: 11 c ------ 12 c 13 c ********************************************************************* 14 c .... calcul du flux de masse aux niveaux s ...... 15 c ********************************************************************* 16 c massebx,masseby,vcont et ucont sont des argum. d'entree pour le s-pg . 17 c pbaru et pbarv sont des argum.de sortie pour le s-pg . 18 c 19 c======================================================================= 37 ijb=ij_begin 38 ije=ij_end+iip1 39 IF(pole_nord) ijb=ij_begin+iip1 40 IF(pole_sud) ije=ij_end-iip1 41 pbaru(ijb:ije,l)=massebx(ijb:ije,l)*ucont(ijb:ije,l) 20 42 43 ijb=ij_begin-iip1 44 ije=ij_end+iip1 45 IF(pole_nord) ijb=ij_begin 46 IF(pole_sud) ije=ij_end-iip1 47 pbarv(ijb:ije,l)=masseby(ijb:ije,l)*vcont(ijb:ije,l) 21 48 22 #include "dimensions.h" 23 #include "paramet.h" 24 #include "comgeom.h" 49 END DO 50 !$OMP END DO NOWAIT 25 51 26 REAL massebx( ijb_u:ije_u,llm ),masseby( ijb_v:ije_v,llm ) , 27 * vcont( ijb_v:ije_v,llm ),ucont( ijb_u:ije_u,llm ), 28 * pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm ) 52 !--- North pole 53 IF(pole_nord) THEN 54 sairen =SUM(aire (1:iim)) 55 saireun=SUM(aireu(1:iim)) 56 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 57 DO l=1,llm 58 ctn=SUM(pbarv(1:iim,l))/sairen 59 pbaru(1,l)= pbarv(1,l)-ctn*aire(1) 60 DO i=2,iim 61 pbaru(i,l)=pbaru(i-1,l)+pbarv(i,l)-ctn*aire(i) 62 END DO 63 apbarun(:)=aireu(1:iim)*pbaru(1:iim,l) 64 ctn0 = -SUM(apbarun)/saireun 65 pbaru(1:iim,l)=2.*(pbaru(1:iim,l)+ctn0) 66 pbaru(iip1,l)=pbaru(1,l) 67 END DO 68 !$OMP END DO NOWAIT 69 END IF 29 70 30 REAL apbarun( iip1 ),apbarus( iip1 ) 71 !--- South pole 72 IF(pole_sud) THEN 73 saires =SUM(aire (ip1jm+1:ip1jmp1-1)) 74 saireus=SUM(aireu(ip1jm+1:ip1jmp1-1)) 75 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 76 DO l=1,llm 77 cts=SUM(pbarv(1+ip1jmi1:ip1jm-1,l))/saires 78 pbaru(1+ip1jm,l)=-pbarv(1+ip1jmi1,l)+cts*aire(1+ip1jm) 79 DO i=2,iim 80 pbaru(i+ip1jm,l)=pbaru(i-1+ip1jm,l)-pbarv(i+ip1jmi1,l)+cts*aire(i+ip1jm) 81 END DO 82 apbarus(:)=aireu(1+ip1jm:ip1jmp1-1)*pbaru(1+ip1jm:ip1jmp1-1,l) 83 cts0 = -SUM(apbarus)/saireus 84 pbaru(1+ip1jm:ip1jmp1-1,l)=2.*(pbaru(1+ip1jm:ip1jmp1-1,l)+cts0) 85 pbaru(ip1jmp1,l)=pbaru(1+ip1jm,l) 86 END DO 87 !$OMP END DO NOWAIT 88 END IF 31 89 32 REAL sairen,saireun,saires,saireus,ctn,cts,ctn0,cts0 33 INTEGER l,ij,i 34 INTEGER ijb,ije 35 36 EXTERNAL SSUM 37 REAL SSUM 38 39 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 40 DO 5 l = 1,llm 90 END SUBROUTINE flumass_loc 41 91 42 ijb=ij_begin43 ije=ij_end+iip144 45 if (pole_nord) ijb=ij_begin+iip146 if (pole_sud) ije=ij_end-iip147 48 DO 1 ij = ijb,ije49 pbaru( ij,l ) = massebx( ij,l ) * ucont( ij,l )50 1 CONTINUE51 52 ijb=ij_begin-iip153 ije=ij_end+iip154 55 if (pole_nord) ijb=ij_begin56 if (pole_sud) ije=ij_end-iip157 58 DO 3 ij = ijb,ije59 pbarv( ij,l ) = masseby( ij,l ) * vcont( ij,l )60 3 CONTINUE61 62 5 CONTINUE63 c$OMP END DO NOWAIT64 c ................................................................65 c calcul de la composante du flux de masse en x aux poles .......66 c ................................................................67 c par la resolution d'1 systeme de 2 equations .68 69 c la premiere equat.decrivant le calcul de la divergence en 1 point i70 c du pole,ce calcul etant itere de i=1 a i=im .71 c c.a.d ,72 c ( ( 0.5*pbaru(i)-0.5*pbaru(i-1) - pbarv(i))/aire(i) =73 c - somme de ( pbarv(n) )/aire pole74 75 c l'autre equat.specifiant que la moyenne du flux de masse au pole est =0.76 c c.a.d somme de pbaru(n)*aire locale(n) = 0.77 78 c on en revient ainsi a determiner la constante additive commune aux pbaru79 c qui representait pbaru(0,j,l) dans l'equat.du calcul de la diverg.au pt80 c i=1 .81 c i variant de 1 a im82 c n variant de 1 a im83 84 IF (pole_nord) THEN85 86 sairen = SSUM( iim, aire( 1 ), 1 )87 saireun= SSUM( iim, aireu( 1 ), 1 )88 89 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)90 DO l = 1,llm91 92 ctn = SSUM( iim, pbarv( 1 ,l), 1 )/ sairen93 94 pbaru(1,l)=pbarv(1,l) - ctn * aire(1)95 96 DO i = 2,iim97 pbaru(i,l) = pbaru(i- 1,l ) +98 * pbarv(i,l) - ctn * aire(i )99 ENDDO100 101 DO i = 1,iim102 apbarun(i) = aireu( i ) * pbaru( i , l)103 ENDDO104 105 ctn0 = -SSUM( iim,apbarun,1 )/saireun106 107 DO i = 1,iim108 pbaru( i , l) = 2. * ( pbaru( i , l) + ctn0 )109 ENDDO110 111 pbaru( iip1 ,l ) = pbaru( 1 ,l )112 113 ENDDO114 c$OMP END DO NOWAIT115 116 ENDIF117 118 119 IF (pole_sud) THEN120 121 saires = SSUM( iim, aire( ip1jm+1 ), 1 )122 saireus= SSUM( iim, aireu( ip1jm+1 ), 1 )123 124 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)125 DO l = 1,llm126 127 cts = SSUM( iim, pbarv(ip1jmi1+ 1,l), 1 )/ saires128 pbaru(ip1jm+1,l)= - pbarv(ip1jmi1+1,l) + cts * aire(ip1jm+1)129 130 DO i = 2,iim131 pbaru(i+ ip1jm,l) = pbaru(i+ip1jm-1,l) -132 * pbarv(i+ip1jmi1,l)+cts*aire(i+ip1jm)133 ENDDO134 135 DO i = 1,iim136 apbarus(i) = aireu(i +ip1jm) * pbaru(i +ip1jm, l)137 ENDDO138 139 cts0 = -SSUM( iim,apbarus,1 )/saireus140 141 DO i = 1,iim142 pbaru(i+ ip1jm, l) = 2. * ( pbaru(i +ip1jm, l) + cts0 )143 ENDDO144 145 pbaru( ip1jmp1,l ) = pbaru( ip1jm +1,l )146 147 ENDDO148 c$OMP END DO NOWAIT149 ENDIF150 151 RETURN152 END -
LMDZ5/trunk/libf/dyn3dmem/massbar_loc.F90
r2335 r2336 1 SUBROUTINE massbar_loc( masse, massebx, masseby ) 2 3 c 4 c ********************************************************************** 5 c 6 c Calcule les moyennes en x et y de la masse d'air dans chaque maille. 7 c ********************************************************************** 8 c Auteurs : P. Le Van , Fr. Hourdin . 9 c .......... 10 c 11 c .. masse est un argum. d'entree pour le s-pg ... 12 c .. massebx,masseby sont des argum. de sortie pour le s-pg ... 13 c 14 c 15 USE parallel_lmdz 16 IMPLICIT NONE 17 c 18 #include "dimensions.h" 19 #include "paramet.h" 20 #include "comconst.h" 21 #include "comgeom.h" 22 c 23 REAL masse( ijb_u:ije_u,llm ), massebx( ijb_u:ije_u,llm ) , 24 * masseby( ijb_v:ije_v,llm ) 25 INTEGER ij,l,ijb,ije 26 c 27 c 28 c Methode pour calculer massebx et masseby . 29 c ---------------------------------------- 30 c 31 c A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires 32 c alpha1(i,j) calcule au point ( i+1/4,j-1/4 ) 33 c alpha2(i,j) calcule au point ( i+1/4,j+1/4 ) 34 c alpha3(i,j) calcule au point ( i-1/4,j+1/4 ) 35 c alpha4(i,j) calcule au point ( i-1/4,j-1/4 ) 36 c 37 c Avec alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j) 38 c 39 c N.B . Pour plus de details, voir s-pg ... iniconst ... 40 c 41 c 42 c 43 c alpha4 . . alpha1 . alpha4 44 c (i,j) (i,j) (i+1,j) 45 c 46 c P . U . . P 47 c (i,j) (i,j) (i+1,j) 48 c 49 c alpha3 . . alpha2 .alpha3 50 c (i,j) (i,j) (i+1,j) 51 c 52 c V . Z . . V 53 c (i,j) 54 c 55 c alpha4 . . alpha1 .alpha4 56 c (i,j+1) (i,j+1) (i+1,j+1) 57 c 58 c P . U . . P 59 c (i,j+1) (i+1,j+1) 60 c 61 c 62 c 63 c On a : 64 c 65 c massebx(i,j) = masse(i ,j) * ( alpha1(i ,j) + alpha2(i,j)) + 66 c masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) ) 67 c localise au point ... U (i,j) ... 68 c 69 c masseby(i,j) = masse(i,j ) * ( alpha2(i,j ) + alpha3(i,j ) + 70 c masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1) 71 c localise au point ... V (i,j) ... 72 c 73 c 74 c======================================================================= 75 76 77 78 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 79 DO 100 l = 1 , llm 80 c 81 ijb=ij_begin 82 ije=ij_end+iip1 83 if (pole_sud) ije=ije-iip1 84 85 DO ij = ijb, ije - 1 86 massebx(ij,l) = masse( ij, l) * alpha1p2( ij ) + 87 * masse(ij+1, l) * alpha3p4(ij+1 ) 88 ENDDO 1 SUBROUTINE massbar_loc(masse,massebx,masseby) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute air mass mean along X and Y in each cell. 7 ! See iniconst for more details. 8 USE parallel_lmdz 9 IMPLICIT NONE 10 include "dimensions.h" 11 include "paramet.h" 12 include "comgeom.h" 13 !=============================================================================== 14 ! Arguments: 15 REAL, INTENT(IN) :: masse (ijb_u:ije_u,llm) 16 REAL, INTENT(OUT) :: massebx(ijb_u:ije_u,llm) 17 REAL, INTENT(OUT) :: masseby(ijb_v:ije_v,llm) 18 !------------------------------------------------------------------------------- 19 ! Method used. Each scalar point is associated to 4 area coefficients: 20 ! * alpha1(i,j) at point ( i+1/4,j-1/4 ) 21 ! * alpha2(i,j) at point ( i+1/4,j+1/4 ) 22 ! * alpha3(i,j) at point ( i-1/4,j+1/4 ) 23 ! * alpha4(i,j) at point ( i-1/4,j-1/4 ) 24 ! where alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j) 25 ! 26 ! alpha4 . . alpha1 . alpha4 27 ! (i,j) (i,j) (i+1,j) 28 ! 29 ! P . U . . P 30 ! (i,j) (i,j) (i+1,j) 31 ! 32 ! alpha3 . . alpha2 .alpha3 33 ! (i,j) (i,j) (i+1,j) 34 ! 35 ! V . Z . . V 36 ! (i,j) 37 ! 38 ! alpha4 . . alpha1 .alpha4 39 ! (i,j+1) (i,j+1) (i+1,j+1) 40 ! 41 ! P . U . . P 42 ! (i,j+1) (i+1,j+1) 43 ! 44 ! 45 ! massebx(i,j) = masse(i ,j) * ( alpha1(i ,j) + alpha2(i,j)) + 46 ! masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) ) 47 ! localized at point ... U (i,j) ... 48 ! 49 ! masseby(i,j) = masse(i,j ) * ( alpha2(i,j ) + alpha3(i,j ) + 50 ! masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1) 51 ! localized at point ... V (i,j) ... 52 !=============================================================================== 53 ! Local variables: 54 INTEGER :: ij, l, ijb, ije 55 !=============================================================================== 56 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 57 DO l=1,llm 58 ijb=ij_begin 59 ije=ij_end+iip1 60 IF(pole_sud) ije=ije-iip1 61 DO ij=ijb,ije-1 62 massebx(ij,l)=masse(ij,l)*alpha1p2(ij)+masse(ij+1 ,l)*alpha3p4(ij+1) 63 END DO 64 DO ij=ijb+iim,ije+iim,iip1; massebx(ij,l)=massebx(ij-iim,l); END DO 65 ijb=ij_begin-iip1 66 ije=ij_end+iip1 67 IF(pole_nord) ijb=ij_begin 68 IF(pole_sud) ije=ij_end-iip1 69 DO ij=ijb,ije 70 masseby(ij,l)=masse(ij,l)*alpha2p3(ij)+masse(ij+iip1,l)*alpha1p4(ij+iip1) 71 END DO 72 END DO 73 !$OMP END DO NOWAIT 89 74 90 c .... correction pour massebx( iip1,j) ..... 91 c ... massebx(iip1,j)= massebx(1,j) ... 92 c 93 CDIR$ IVDEP 75 END SUBROUTINE massbar_loc 94 76 95 96 97 DO ij = ijb+iim, ije+iim, iip198 massebx( ij,l ) = massebx( ij - iim,l )99 ENDDO100 101 102 103 ijb=ij_begin-iip1104 ije=ij_end+iip1105 if (pole_nord) ijb=ij_begin106 if (pole_sud) ije=ij_end-iip1107 108 DO ij = ijb,ije109 masseby( ij,l ) = masse( ij , l ) * alpha2p3( ij ) +110 * masse(ij+iip1, l ) * alpha1p4( ij+iip1 )111 ENDDO112 113 100 CONTINUE114 c$OMP END DO NOWAIT115 c116 RETURN117 END -
LMDZ5/trunk/libf/dyn3dmem/massbarxy_loc.F90
r2335 r2336 1 SUBROUTINE massbarxy_loc( masse, massebxy ) 2 USE parallel_lmdz 3 implicit none 4 c ********************************************************************** 5 c 6 c Calcule les moyennes en x et y de la masse d'air dans chaque maille. 7 c ********************************************************************** 8 c Auteurs : P. Le Van , Fr. Hourdin . 9 c .......... 10 c 11 c .. masse est un argum. d'entree pour le s-pg ... 12 c .. massebxy est un argum. de sortie pour le s-pg ... 13 c 14 c 15 c IMPLICIT NONE 16 c 17 #include "dimensions.h" 18 #include "paramet.h" 19 #include "comconst.h" 20 #include "comgeom.h" 21 c 22 REAL masse( ijb_u:ije_u,llm ), massebxy( ijb_v:ije_v,llm ) 23 c 24 INTEGER ij,l,ijb,ije 1 SUBROUTINE massbarxy_loc(masse,massebxy) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute air mass mean along X and Y in each cell. 7 ! See iniconst for more details. 8 USE parallel_lmdz 9 IMPLICIT NONE 10 include "dimensions.h" 11 include "paramet.h" 12 include "comconst.h" 13 include "comgeom.h" 14 !=============================================================================== 15 ! Arguments: 16 REAL, INTENT(IN) :: masse (ijb_u:ije_u,llm) 17 REAL, INTENT(OUT) :: massebxy(ijb_v:ije_v,llm) 18 !=============================================================================== 19 ! Local variables: 20 INTEGER :: ij, l, ijb, ije 21 !=============================================================================== 22 ijb=ij_begin-iip1 23 ije=ij_end 24 IF(pole_nord) ijb=ijb+iip1 25 IF(pole_sud) ije=ije-iip1 26 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 27 DO l=1,llm 28 DO ij=ijb,ije-1 29 massebxy(ij,l)=masse(ij ,l)*alpha2(ij ) + & 30 + masse(ij+1 ,l)*alpha3(ij+1 ) + & 31 + masse(ij+iip1,l)*alpha1(ij+iip1) + & 32 + masse(ij+iip2,l)*alpha4(ij+iip2) 33 END DO 34 DO ij=ijb+iip1-1,ije+iip1-1,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO 35 END DO 36 !$OMP END DO NOWAIT 25 37 26 27 ijb=ij_begin-iip1 28 ije=ij_end 29 30 if (pole_nord) ijb=ijb+iip1 31 if (pole_sud) ije=ije-iip1 38 END SUBROUTINE massbarxy_loc 32 39 33 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)34 DO 100 l = 1 , llm35 c36 DO 5 ij = ijb, ije - 137 massebxy( ij,l ) = masse( ij ,l ) * alpha2( ij ) +38 + masse( ij+1 ,l ) * alpha3( ij+1 ) +39 + masse( ij+iip1,l ) * alpha1( ij+iip1 ) +40 + masse( ij+iip2,l ) * alpha4( ij+iip2 )41 5 CONTINUE42 43 c .... correction pour massebxy( iip1,j ) ........44 45 CDIR$ IVDEP46 47 DO 7 ij = ijb+iip1-1, ije+iip1-1, iip148 massebxy( ij,l ) = massebxy( ij - iim,l )49 7 CONTINUE50 51 100 CONTINUE52 c$OMP END DO NOWAIT53 c54 RETURN55 END -
LMDZ5/trunk/libf/dyn3dmem/tourpot_loc.F90
r2335 r2336 1 SUBROUTINE tourpot_loc ( vcov, ucov, massebxy, vorpot ) 2 USE parallel_lmdz 3 USE mod_filtreg_p 4 IMPLICIT NONE 1 SUBROUTINE tourpot_loc ( vcov, ucov, massebxy, vorpot ) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute potential vorticity. 7 USE parallel_lmdz 8 USE mod_filtreg_p 9 IMPLICIT NONE 10 include "dimensions.h" 11 include "paramet.h" 12 include "comgeom.h" 13 include "logic.h" 14 !=============================================================================== 15 ! Arguments: 16 REAL, INTENT(IN) :: vcov (ijb_v:ije_v,llm) 17 REAL, INTENT(IN) :: ucov (ijb_u:ije_u,llm) 18 REAL, INTENT(IN) :: massebxy(ijb_v:ije_v,llm) 19 REAL, INTENT(OUT) :: vorpot (ijb_v:ije_v,llm) 20 !=============================================================================== 21 ! Method used: 22 ! vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy 23 !=============================================================================== 24 ! Local variables: 25 INTEGER :: l, ij, ije, ijb, jje, jjb 26 REAL :: rot(ijb_v:ije_v,llm) 27 !=============================================================================== 5 28 6 c======================================================================= 7 c 8 c Auteur: P. Le Van 9 c ------- 10 c 11 c Objet: 12 c ------ 13 c 14 c ******************************************************************* 15 c ......... calcul du tourbillon potentiel ......... 16 c ******************************************************************* 17 c 18 c vcov,ucov,fext et pbarxyfl sont des argum. d'entree pour le s-pg . 19 c vorpot est un argum.de sortie pour le s-pg . 20 c 21 c======================================================================= 29 ijb=ij_begin-iip1 30 ije=ij_end 31 IF(pole_nord) ijb=ij_begin 22 32 23 #include "dimensions.h" 24 #include "paramet.h" 25 #include "comgeom.h" 26 #include "logic.h" 33 !--- Wind vorticity ; correction: rot(iip1,j,l) = rot(1,j,l) 34 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 35 DO l=1,llm 36 IF(pole_sud) ije=ij_end-iip1-1 37 DO ij=ijb,ije 38 rot(ij,l)=vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l) 39 END DO 40 IF(pole_sud) ije=ij_end-iip1 41 DO ij=ijb+iip1-1,ije,iip1; rot(ij,l)=rot(ij-iim,l); END DO 42 END DO 43 !$OMP END DO NOWAIT 27 44 28 REAL rot( ijb_v:ije_v,llm ) 29 REAL vcov( ijb_v:ije_v,llm ),ucov( ijb_u:ije_u,llm ) 30 REAL massebxy( ijb_v:ije_v,llm ),vorpot( ijb_v:ije_v,llm ) 45 !--- Filter 46 jjb=jj_begin-1 47 jje=jj_end 48 IF(pole_nord) jjb=jjb+1 49 IF(pole_sud) jje=jje-1 50 CALL filtreg_p(rot,jjb_v,jje_v,jjb,jje,jjm,llm,2,1,.FALSE.,1) 31 51 32 INTEGER l, ij ,ije,ijb,jje,jjb 52 !--- Potential vorticity ; correction: rot(iip1,j,l) = rot(1,j,l) 53 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 54 DO l=1,llm 55 IF(pole_sud) ije=ij_end-iip1-1 56 DO ij=ijb,ije 57 vorpot(ij,l)=(rot(ij,l)+fext(ij))/massebxy(ij,l) 58 END DO 59 IF(pole_sud) ije=ij_end-iip1 60 DO ij=ijb+iip1-1,ije,iip1; vorpot(ij,l)=vorpot(ij-iim,l); END DO 61 END DO 62 !$OMP END DO NOWAIT 33 63 64 END SUBROUTINE tourpot_loc 34 65 35 ijb=ij_begin-iip136 ije=ij_end37 38 if (pole_nord) ijb=ij_begin39 40 41 c ... vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy ..42 43 44 45 c ........ Calcul du rotationnel du vent V puis filtrage ........46 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)47 DO 5 l = 1,llm48 49 if (pole_sud) ije=ij_end-iip1-150 DO 2 ij = ijb, ije51 rot( ij,l ) = vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)52 2 CONTINUE53 54 c .... correction pour rot( iip1,j,l ) .....55 c .... rot(iip1,j,l) = rot(1,j,l) .....56 57 CDIR$ IVDEP58 59 if (pole_sud) ije=ij_end-iip160 61 DO 3 ij = ijb+iip1-1, ije, iip162 rot( ij,l ) = rot( ij -iim, l )63 3 CONTINUE64 65 5 CONTINUE66 c$OMP END DO NOWAIT67 jjb=jj_begin-168 jje=jj_end69 70 if (pole_nord) jjb=jjb+171 if (pole_sud) jje=jje-172 CALL filtreg_p( rot, jjb_v,jje_v,jjb,jje,jjm, llm,73 & 2, 1, .FALSE., 1 )74 75 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)76 DO 10 l = 1, llm77 78 if (pole_sud) ije=ij_end-iip1-179 80 DO 6 ij = ijb, ije81 vorpot( ij,l ) = ( rot(ij,l) + fext(ij) ) / massebxy(ij,l)82 6 CONTINUE83 84 c ..... correction pour vorpot( iip1,j,l) .....85 c .... vorpot(iip1,j,l)= vorpot(1,j,l) ....86 CDIR$ IVDEP87 if (pole_sud) ije=ij_end-iip188 DO 8 ij = ijb+iip1-1, ije, iip189 vorpot( ij,l ) = vorpot( ij -iim,l )90 8 CONTINUE91 92 10 CONTINUE93 c$OMP END DO NOWAIT94 RETURN95 END -
LMDZ5/trunk/libf/dyn3dmem/vitvert_loc.F90
r2335 r2336 1 SUBROUTINE vitvert_loc ( convm , w ) 2 c 3 USE parallel_lmdz 4 IMPLICIT NONE 1 SUBROUTINE vitvert_loc(convm, w) 2 ! 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr. Hourdin. 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute vertical speed at sigma levels. 7 USE parallel_lmdz 8 IMPLICIT NONE 9 include "dimensions.h" 10 include "paramet.h" 11 include "comvert.h" 12 !=============================================================================== 13 ! Arguments: 14 REAL, INTENT(IN) :: convm(ijb_u:ije_u,llm) 15 REAL, INTENT(OUT) :: w (ijb_u:ije_u,llm) 16 !=============================================================================== 17 ! Notes: Vertical speed is oriented from bottom to top. 18 ! * At ground - level sigma(1): w(i,j,1) = 0. 19 ! * At top - level sigma(llm+1): w(i,j,l) = 0. (not stored in w) 20 !=============================================================================== 21 ! Local variables: 22 INTEGER :: l, ijb, ije 23 !=============================================================================== 24 ijb=ij_begin 25 ije=ij_end+iip1 26 IF(pole_sud) ije=ij_end 27 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 28 DO l=1,llmm1 29 w(ijb:ije,l+1)=convm(ijb:ije,l+1)-bp(l+1)*convm(ijb:ije,1) 30 END DO 31 !$OMP END DO 32 !$OMP MASTER 33 w(ijb:ije,1)=0. 34 !$OMP END MASTER 35 !$OMP BARRIER 5 36 6 c======================================================================= 7 c 8 c Auteurs: P. Le Van , F. Hourdin . 9 c ------- 10 c 11 c Objet: 12 c ------ 13 c 14 c ******************************************************************* 15 c .... calcul de la vitesse verticale aux niveaux sigma .... 16 c ******************************************************************* 17 c convm est un argument d'entree pour le s-pg ...... 18 c w est un argument de sortie pour le s-pg ...... 19 c 20 c la vitesse verticale est orientee de haut en bas . 21 c au sol, au niveau sigma(1), w(i,j,1) = 0. 22 c au sommet, au niveau sigma(llm+1) , la vit.verticale est aussi 23 c egale a 0. et n'est pas stockee dans le tableau w . 24 c 25 c 26 c======================================================================= 37 END SUBROUTINE vitvert_loc 27 38 28 #include "dimensions.h"29 #include "paramet.h"30 #include "comvert.h"31 32 REAL w(ijb_u:ije_u,llm),convm(ijb_u:ije_u,llm)33 INTEGER l, ij,ijb,ije34 35 36 ijb=ij_begin37 ije=ij_end+iip138 39 if (pole_sud) ije=ij_end40 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)41 DO 2 l = 1,llmm142 43 DO 1 ij = ijb,ije44 w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )45 1 CONTINUE46 47 2 CONTINUE48 c$OMP END DO49 c$OMP MASTER50 DO 5 ij = ijb,ije51 w(ij,1) = 0.52 5 CONTINUE53 c$OMP END MASTER54 c$OMP BARRIER55 RETURN56 END -
LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90
r2331 r2336 1 1 PROGRAM ce0l 2 ! 3 ! Purpose: Calls etat0, creates initial states and limit_netcdf 4 ! 5 ! interbar=.T. for barycentric interpolation inter_barxy 6 ! extrap =.T. for data extrapolation, like for the SSTs when file does not 7 ! contain ocean points only. 8 ! oldice =.T. for old-style ice, obtained using grille_m (grid_atob). 9 ! masque is created in etat0, passed to limit to ensure consistancy. 10 11 USE control_mod, only: DAY_STEP, DAYREF, NSPLIT_PHYS 12 USE etat0dyn, only: etat0dyn_netcdf 13 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR 2 ! 3 !------------------------------------------------------------------------------- 4 ! Purpose: Initial states and boundary conditions files creation: 5 ! * start.nc for dynamics (using etat0dyn routine) 6 ! * startphy.nc for physics (using etat0phys routine) 7 ! * limit.nc for forced runs (using limit_netcdf routine) 8 !------------------------------------------------------------------------------- 9 ! Notes: 10 ! * extrap=.T. (default) for data extrapolation, like for the SSTs when file 11 ! does contain ocean points only. 12 ! * "masque" can be: 13 ! - read from file "o2a.nc" (for coupled runs). 14 ! - created in etat0phys or etat0dyn (for forced runs). 15 ! It is then passed to limit_netcdf to ensure consistancy. 16 !------------------------------------------------------------------------------- 14 17 USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo 15 16 USE etat0phys, only: etat0phys_netcdf 17 USE dimphy, only: KLON 18 USE infotrac, only: TYPE_TRAC, infotrac_init 18 USE control_mod, ONLY: day_step, dayref, nsplit_phys 19 USE etat0dyn, ONLY: etat0dyn_netcdf 20 USE etat0phys, ONLY: etat0phys_netcdf 21 USE limit, ONLY: limit_netcdf 22 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR 23 USE infotrac, ONLY: type_trac, infotrac_init 24 USE dimphy, ONLY: klon 19 25 USE test_disvert_m, ONLY: test_disvert 26 USE filtreg_mod, ONLY: inifilr 27 #ifdef inca 28 USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic 29 #endif 30 #ifdef CPP_PARA 31 USE mod_const_mpi, ONLY: init_const_mpi 32 USE parallel_lmdz, ONLY: init_parallel, mpi_rank, omp_rank, mpi_size 33 USE bands, ONLY: read_distrib, distrib_phys 34 USE mod_hallo, ONLY: init_mod_hallo 35 USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys 36 #endif 20 37 21 38 IMPLICIT NONE 22 39 23 ! Local variables: 40 !------------------------------------------------------------------------------- 41 ! Local variables: 24 42 include "dimensions.h" 25 43 include "paramet.h" 26 include "comgeom .h"44 include "comgeom2.h" 27 45 include "comconst.h" 28 46 include "comvert.h" … … 30 48 include "temps.h" 31 49 include "logic.h" 32 REAL :: masque(iip1, 33 REAL :: phis (iip1, 50 REAL :: masque(iip1,jjp1) !--- CONTINENTAL MASK 51 REAL :: phis (iip1,jjp1) !--- GROUND GEOPOTENTIAL 34 52 CHARACTER(LEN=256) :: modname, fmt, calnd !--- CALENDAR TYPE 35 53 LOGICAL :: use_filtre_fft 36 LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.37 38 54 LOGICAL, PARAMETER :: extrap=.FALSE. 55 56 !--- Local variables for ocean mask reading: 39 57 INTEGER :: nid_o2a, iml_omask, jml_omask, j 40 58 INTEGER :: fid, iret, llm_tmp, ttm_tmp, itaul(1) 41 REAL, ALLOCATABLE :: lon_omask(:, :), dlon_omask(:), ocemask(:,:)42 REAL, ALLOCATABLE :: lat_omask(:, :), dlat_omask(:), ocetmp (:,:)59 REAL, ALLOCATABLE :: lon_omask(:,:), dlon_omask(:), ocemask(:,:) 60 REAL, ALLOCATABLE :: lat_omask(:,:), dlat_omask(:), ocetmp (:,:) 43 61 REAL :: date, lev(1) 44 45 !---------------------------------------------------------------------- 62 !------------------------------------------------------------------------------- 46 63 modname="ce0l" 47 64 48 65 !--- Constants 49 66 pi = 4. * ATAN(1.) 50 67 rad = 6371229. … … 59 76 60 77 CALL conf_gcm( 99, .TRUE. ) 61 62 78 dtvr = daysec/REAL(day_step) 63 WRITE(lunout, *)'dtvr', dtvr 64 79 WRITE(lunout,*)'dtvr',dtvr 65 80 CALL iniconst() 66 81 CALL inigeom() 67 82 83 !--- Calendar choice 68 84 #ifdef CPP_IOIPSL 69 85 calnd='gregorian' 70 86 SELECT CASE(calend) 71 CASE('earth_360d') 72 CALL ioconf_calendar('360d') 73 calnd='with 360 days/year' 74 CASE('earth_365d') 75 CALL ioconf_calendar('noleap') 76 calnd='with no leap year' 77 CASE('earth_366d') 78 CALL ioconf_calendar('366d') 79 calnd='with leap years only' 80 CASE('gregorian') 81 CALL ioconf_calendar('gregorian') 82 CASE('standard') 83 CALL ioconf_calendar('gregorian') 84 CASE('julian') 85 CALL ioconf_calendar('julian') 86 calnd='julian' 87 CASE('proleptic_gregorian') 88 CALL ioconf_calendar('gregorian') 89 !--- DC Bof... => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian 90 CASE DEFAULT 91 CALL abort_gcm('ce0l', 'Bad choice for calendar', 1) 87 CASE('earth_360d');CALL ioconf_calendar('360d'); calnd='with 360 days/year' 88 CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='with no leap year' 89 CASE('earth_366d');CALL ioconf_calendar('366d'); calnd='with leap years only' 90 CASE('gregorian'); CALL ioconf_calendar('gregorian') 91 CASE('standard'); CALL ioconf_calendar('gregorian') 92 CASE('julian'); CALL ioconf_calendar('julian'); calnd='julian' 93 CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian') 94 !--- DC Bof... => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian 95 CASE DEFAULT 96 CALL abort_gcm('ce0l','Bad choice for calendar',1) 92 97 END SELECT 93 WRITE(lunout, *)'CHOSEN CALENDAR: Earth '//TRIM(calnd) 94 #endif 95 96 use_filtre_fft=.FALSE. 97 CALL getin('use_filtre_fft', use_filtre_fft) 98 IF(use_filtre_fft) THEN 99 WRITE(lunout, *)"FFT filter not available for sequential dynamics." 100 WRITE(lunout, *)"Your setting of variable use_filtre_fft is not used." 101 ENDIF 102 103 !--- LAND MASK. TWO CASES: 104 ! 1) read from ocean model file "o2a.nc" (coupled runs) 105 ! 2) computed from topography file "Relief.nc" (masque(:, :)=-99999.) 106 ! Coupled simulations (case 1) use the ocean model mask to compute the 107 ! weights to ensure ocean fractions are the same for atmosphere and ocean. 108 109 IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN 110 WRITE(lunout, *)'BEWARE !! No ocean mask "o2a.nc" file found' 111 WRITE(lunout, *)'Forced run.' 112 masque(:, :)=-99999. 113 ELSE 114 iret=NF90_CLOSE(nid_o2a) 115 WRITE(lunout, *)'BEWARE !! Ocean mask "o2a.nc" file found' 116 WRITE(lunout, *)'Coupled run.' 117 CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a) 118 IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN 119 WRITE(lunout, *)'Mismatching dimensions for ocean mask' 120 WRITE(lunout, *)'iim = ', iim , ' iml_omask = ', iml_omask 121 WRITE(lunout, *)'jjp1 = ', jjp1, ' jml_omask = ', jml_omask 122 CALL abort_gcm(modname, '', 1) 123 END IF 124 ALLOCATE(ocemask(iim, jjp1), lon_omask(iim, jjp1), dlon_omask(iim )) 125 ALLOCATE(ocetmp (iim, jjp1), lat_omask(iim, jjp1), dlat_omask(jjp1)) 126 CALL flinopen("o2a.nc", .FALSE., iim, jjp1, llm_tmp, lon_omask, & 127 lat_omask, lev, ttm_tmp, itaul, date, dt, fid) 128 CALL flinget(fid, "OceMask", iim, jjp1, llm_tmp, ttm_tmp, 1, 1, ocetmp) 129 CALL flinclo(fid) 130 dlon_omask(1:iim ) = lon_omask(1:iim, 1) 131 dlat_omask(1:jjp1) = lat_omask(1, 1:jjp1) 132 ocemask = ocetmp 133 IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN 134 DO j=1, jjp1 135 ocemask(:, j) = ocetmp(:, jjp1-j+1) 136 END DO 137 END IF 138 DEALLOCATE(ocetmp, lon_omask, lat_omask, dlon_omask, dlat_omask) 139 IF(prt_level>=1) THEN 140 WRITE(fmt, "(i4, 'i1)')")iim 141 fmt='('//ADJUSTL(fmt) 142 WRITE(lunout, *)'OCEAN MASK :' 143 WRITE(lunout, fmt) NINT(ocemask) 144 END IF 145 masque(1:iim, :)=1.-ocemask(:, :) 146 masque(iip1 , :)=masque(1, :) 147 DEALLOCATE(ocemask) 148 END IF 149 phis(:, :)=-99999. 150 151 CALL Init_Phys_lmdz(iim, jjp1, llm, 1, (/(jjm-1)*iim+2/)) 152 WRITE(lunout, *)'---> klon=', klon 153 154 call infotrac_init 155 CALL iniphysiq(iim, jjm, llm, daysec, dayref, dtphys / nsplit_phys, rlatu, & 156 rlonv, aire, cu, cv, rad, g, r, cpp, iflag_phys) 157 158 IF(pressure_exner) CALL test_disvert 159 98 WRITE(lunout,*)'CHOSEN CALENDAR: Earth '//TRIM(calnd) 99 #endif 100 101 #ifdef CPP_PARA 102 !--- Physical grid + parallel initializations 103 CALL init_const_mpi() 104 CALL init_parallel() 105 CALL read_distrib() 106 CALL init_mod_hallo() 107 CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys) 108 CALL init_interface_dyn_phys() 109 #else 110 CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) 111 #endif 112 WRITE(lunout,*)'---> klon=',klon 113 114 !--- Tracers initializations 160 115 IF (type_trac == 'inca') THEN 161 116 #ifdef INCA 162 CALL init_const_lmdz(nbtr, anneeref, dayref, iphysiq, day_step, nday) 163 CALL init_inca_para(iim, jjm+1, klon, 1, klon_mpi_para_nb, 0) 164 WRITE(lunout, *)'nbtr =' , nbtr 165 #endif 166 END IF 117 CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,& 118 nbsrf,is_oce,is_sic,is_ter,is_lic,calend) 119 CALL init_inca_para(iim,jjp1,llm,klon_glo,mpi_size,distrib_phys,& 120 COMM_LMDZ) 121 WRITE(lunout,*)'nbtr =' , nbtr 122 #endif 123 END IF 124 CALL infotrac_init() 125 126 CALL inifilr() 127 CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, & 128 rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys) 129 IF(pressure_exner) CALL test_disvert 130 131 #ifdef CPP_PARA 132 IF (mpi_rank==0.AND.omp_rank==0) THEN 133 #endif 134 use_filtre_fft=.FALSE. 135 CALL getin('use_filtre_fft',use_filtre_fft) 136 IF(use_filtre_fft) THEN 137 WRITE(lunout,*)"FFT filter not available for sequential dynamics." 138 WRITE(lunout,*)"Your setting of variable use_filtre_fft is not used." 139 ENDIF 140 141 !--- LAND MASK. TWO CASES: 142 ! 1) read from ocean model file "o2a.nc" (coupled runs) 143 ! 2) computed from topography file "Relief.nc" (masque(:,:)=-99999.) 144 ! Coupled simulations (case 1) use the ocean model mask to compute the 145 ! weights to ensure ocean fractions are the same for atmosphere and ocean. 146 !******************************************************************************* 147 IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN 148 WRITE(lunout,*)'BEWARE !! No ocean mask "o2a.nc" file found' 149 WRITE(lunout,*)'Forced run.' 150 masque(:,:)=-99999. 151 ELSE 152 iret=NF90_CLOSE(nid_o2a) 153 WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found' 154 WRITE(lunout,*)'Coupled run.' 155 CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a) 156 IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN 157 WRITE(lunout,*)'Mismatching dimensions for ocean mask' 158 WRITE(lunout,*)'iim = ',iim ,' iml_omask = ',iml_omask 159 WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask 160 CALL abort_gcm(modname,'',1) 161 END IF 162 ALLOCATE(ocemask(iim,jjp1),lon_omask(iim,jjp1),dlon_omask(iim )) 163 ALLOCATE(ocetmp (iim,jjp1),lat_omask(iim,jjp1),dlat_omask(jjp1)) 164 CALL flinopen("o2a.nc", .FALSE.,iim,jjp1,llm_tmp,lon_omask,lat_omask,lev, & 165 & ttm_tmp,itaul,date,dt,fid) 166 CALL flinget(fid, "OceMask", iim,jjp1,llm_tmp,ttm_tmp,1,1,ocetmp) 167 CALL flinclo(fid) 168 dlon_omask(1:iim ) = lon_omask(1:iim,1) 169 dlat_omask(1:jjp1) = lat_omask(1,1:jjp1) 170 ocemask = ocetmp 171 IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN 172 DO j=1,jjp1; ocemask(:,j) = ocetmp(:,jjp1-j+1); END DO 173 END IF 174 DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask) 175 IF(prt_level>=1) THEN 176 WRITE(fmt,"(i4,'i1)')")iim ; fmt='('//ADJUSTL(fmt) 177 WRITE(lunout,*)'OCEAN MASK :' 178 WRITE(lunout,fmt) NINT(ocemask) 179 END IF 180 masque(1:iim,:)=1.-ocemask(:,:) 181 masque(iip1 ,:)=masque(1,:) 182 DEALLOCATE(ocemask) 183 END IF 184 phis(:,:)=-99999. 185 167 186 IF(ok_etat0) THEN 168 WRITE(lunout, '(//)') 169 WRITE(lunout, *) ' ************************ ' 170 WRITE(lunout, *) ' *** etat0phy_netcdf *** ' 171 WRITE(lunout, *) ' ************************ ' 172 WRITE(lunout, '(//)') 173 WRITE(lunout, *) ' interbar = ', interbar 174 CALL etat0phys_netcdf(interbar, masque, phis) 175 END IF 176 177 IF(ok_etat0) THEN 178 WRITE(lunout, '(//)') 179 WRITE(lunout, *) ' ************************ ' 180 WRITE(lunout, *) ' *** etat0dyn_netcdf *** ' 181 WRITE(lunout, *) ' ************************ ' 182 WRITE(lunout, '(//)') 183 WRITE(lunout, *) ' interbar = ', interbar 184 CALL etat0dyn_netcdf(interbar, masque, phis) 187 WRITE(lunout,'(//)') 188 WRITE(lunout,*) ' ************************ ' 189 WRITE(lunout,*) ' *** etat0phy_netcdf *** ' 190 WRITE(lunout,*) ' ************************ ' 191 CALL etat0phys_netcdf(masque,phis) 192 WRITE(lunout,'(//)') 193 WRITE(lunout,*) ' ************************ ' 194 WRITE(lunout,*) ' *** etat0dyn_netcdf *** ' 195 WRITE(lunout,*) ' ************************ ' 196 CALL etat0dyn_netcdf(masque,phis) 185 197 END IF 186 198 187 199 IF(ok_limit) THEN 188 WRITE(lunout, '(//)') 189 WRITE(lunout, *) ' ********************* ' 190 WRITE(lunout, *) ' *** Limit_netcdf *** ' 191 WRITE(lunout, *) ' ********************* ' 192 WRITE(lunout, '(//)') 193 CALL limit_netcdf(interbar, extrap, oldice, masque) 194 END IF 195 196 WRITE(lunout, '(//)') 197 WRITE(lunout, *) ' *************************** ' 198 WRITE(lunout, *) ' *** grilles_gcm_netcdf *** ' 199 WRITE(lunout, *) ' *************************** ' 200 WRITE(lunout, '(//)') 201 CALL grilles_gcm_netcdf_sub(masque, phis) 200 WRITE(lunout,'(//)') 201 WRITE(lunout,*) ' ********************* ' 202 WRITE(lunout,*) ' *** Limit_netcdf *** ' 203 WRITE(lunout,*) ' ********************* ' 204 WRITE(lunout,'(//)') 205 CALL limit_netcdf(masque,phis,extrap) 206 END IF 207 208 WRITE(lunout,'(//)') 209 WRITE(lunout,*) ' *************************** ' 210 WRITE(lunout,*) ' *** grilles_gcm_netcdf *** ' 211 WRITE(lunout,*) ' *************************** ' 212 WRITE(lunout,'(//)') 213 CALL grilles_gcm_netcdf_sub(masque,phis) 214 215 #ifdef CPP_PARA 216 END IF 217 #endif 202 218 203 219 END PROGRAM ce0l 220 ! 221 !------------------------------------------------------------------------------- -
LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90
r2334 r2336 12 12 ! routine (to be called after restget): 13 13 ! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,& 14 ! champ, lon_in2, lat_in2 , ibar)14 ! champ, lon_in2, lat_in2) 15 15 ! 16 16 ! * Variables should have the following names in the NetCDF files: … … 36 36 USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, histclo 37 37 USE assert_eq_m, ONLY: assert_eq 38 #ifdef CPP_PHYS39 38 USE indice_sol_mod, ONLY: epsfra 40 #endif41 39 IMPLICIT NONE 42 40 … … 54 52 include "serre.h" 55 53 REAL, SAVE :: deg2rad 56 #ifndef CPP_PHYS57 REAL, SAVE :: epsfra= 1.E-558 #endif59 54 INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn 60 55 REAL, ALLOCATABLE, SAVE :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:) 61 56 CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc' 62 CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc'63 57 64 58 CONTAINS … … 66 60 !------------------------------------------------------------------------------- 67 61 ! 68 SUBROUTINE etat0dyn_netcdf( ib,masque, phis)62 SUBROUTINE etat0dyn_netcdf(masque, phis) 69 63 ! 70 64 !------------------------------------------------------------------------------- … … 73 67 ! Notes: 1) This routine is designed to work for Earth 74 68 ! 2) If masque(:,:)/=-99999., masque and phis are already known. 75 ! Otherwise: read it from ocean model file (coupled run) orcompute it.69 ! Otherwise: compute it. 76 70 !------------------------------------------------------------------------------- 77 71 USE control_mod 78 #ifdef CPP_PHYS79 72 USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz 80 73 USE regr_pr_o3_m, ONLY: regr_pr_o3 81 74 USE press_coefoz_m, ONLY: press_coefoz 82 #endif83 75 USE exner_hyb_m, ONLY: exner_hyb 84 76 USE exner_milieu_m, ONLY: exner_milieu 85 USE infotrac, only: NQTOT, TNAME77 USE infotrac, ONLY: nqtot, tname 86 78 USE filtreg_mod 87 79 IMPLICIT NONE 88 80 !------------------------------------------------------------------------------- 89 81 ! Arguments: 90 LOGICAL, INTENT(IN) :: ib !--- Barycentric interpolation91 82 REAL, INTENT(INOUT) :: masque(iip1,jjp1) !--- Land-ocean mask 92 83 REAL, INTENT(INOUT) :: phis (iip1,jjp1) !--- Ground geopotential … … 111 102 deg2rad = pi/180.0 112 103 113 ! Initializations for tracers and filter114 !*******************************************************************************115 CALL inifilr()116 117 104 ! Compute ground geopotential and possibly the mask. 118 105 !******************************************************************************* 119 106 masque_tmp(:,:)=masque(:,:) 120 CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp)121 107 WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt) 122 108 IF(ALL(masque==-99999.)) THEN !--- KEEP NEW MASK … … 132 118 ! Compute psol AND tsol, knowing phis. 133 119 !******************************************************************************* 134 CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, ib,phis, psol)120 CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol) 135 121 136 122 ! Mid-levels pressure computation … … 147 133 !******************************************************************************* 148 134 uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0 149 CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv ,ib)135 CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv) 150 136 CALL startget_dyn3d('v' ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent, & 151 & rlonu,rlatu(:jjm) ,ib)152 CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv ,ib)137 & rlonu,rlatu(:jjm)) 138 CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv) 153 139 tpot(:,:,:)=t3d(:,:,:) 154 CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv ,ib)140 CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv) 155 141 156 142 WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:)) … … 165 151 ! WRITE(lunout,*) 'QSAT :',qsat(10,20,:) 166 152 qd (:,:,:) = 0.0 167 CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv ,ib)153 CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv) 168 154 ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:) 169 155 CALL flinclo(fid_dyn) 170 156 171 #ifdef CPP_PHYS172 157 ! Parameterization of ozone chemistry: 173 158 !******************************************************************************* … … 180 165 q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29. !--- Mole->mass fraction 181 166 END IF 182 183 #endif184 167 q3d(iip1,:,:,:)=q3d(1,:,:,:) 185 186 ! Intermediate computation187 !*******************************************************************************188 CALL massdair(p3d,masse)189 WRITE(lunout,*)' ALPHAX ',alphax190 DO l=1,llm191 xppn(:)=aire(1:iim,1 )*masse(1:iim,1 ,l)192 xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)193 xpn=SUM(xppn)/apoln194 xps=SUM(xpps)/apols195 masse(:,1 ,l)=xpn196 masse(:,jjp1,l)=xps197 END DO198 168 199 169 ! Writing … … 215 185 CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi ) 216 186 WRITE(lunout,*)'sortie geopot' 217 CALL caldyn0( itau, uvent, vvent, tpot, psol, pk, phis, &218 phi, w, time+iday-dayref)187 CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, & 188 phi, w, pbaru, pbarv, time+iday-dayref) 219 189 WRITE(lunout,*)'sortie caldyn0' 190 #ifdef CPP_PARA 191 CALL dynredem0_loc( "start.nc", dayref, phis) 192 #else 220 193 CALL dynredem0( "start.nc", dayref, phis) 194 #endif 221 195 WRITE(lunout,*)'sortie dynredem0' 196 #ifdef CPP_PARA 197 CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) 198 #else 222 199 CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) 200 #endif 223 201 WRITE(lunout,*)'sortie dynredem1' 224 202 CALL histclo() … … 232 210 ! 233 211 SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar,& 234 champ, lon_in2, lat_in2 , ibar)212 champ, lon_in2, lat_in2) 235 213 !------------------------------------------------------------------------------- 236 214 IMPLICIT NONE … … 252 230 REAL, INTENT(IN) :: lon_in2(:) ! dim (iml) 253 231 REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2) 254 LOGICAL, INTENT(IN) :: ibar255 232 !------------------------------------------------------------------------------- 256 233 ! Local variables: … … 287 264 !--- INTERPOLATE 3D FIELD IF NEEDED 288 265 IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2, & 289 lat_in2,pls,champ ,ibar)266 lat_in2,pls,champ) 290 267 291 268 !--- COMPUTE THE REQUIRED FILED … … 317 294 !------------------------------------------------------------------------------- 318 295 ! 319 SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque) 320 ! 321 !------------------------------------------------------------------------------- 322 USE conf_dat_m, ONLY: conf_dat2d 323 IMPLICIT NONE 324 !=============================================================================== 325 ! Purpose: Compute "phis" just like it would be in start_init_orog. 326 !=============================================================================== 327 ! Arguments: 328 REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) 329 REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml) 330 !------------------------------------------------------------------------------- 331 ! Local variables: 332 CHARACTER(LEN=256) :: modname="start_init_orog0" 333 CHARACTER(LEN=256) :: title="RELIEF" 334 INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1) 335 REAL :: lev(1), date, dt 336 REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:) 337 REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:) 338 !------------------------------------------------------------------------------- 339 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml") 340 jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml") 341 IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1) 342 IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1) 343 pi=2.0*ASIN(1.0); deg2rad=pi/180.0 344 IF(ANY(phis/=-99999.)) RETURN !--- phis ALREADY KNOWN 345 346 !--- HIGH RESOLUTION OROGRAPHY 347 CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) 348 349 ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) 350 CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,& 351 lev, ttm_tmp, itau, date, dt, fid) 352 ALLOCATE(relief_hi(iml_rel,jml_rel)) 353 CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) 354 CALL flinclo(fid) 355 356 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 357 ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) 358 lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad 359 lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad 360 361 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 362 ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) 363 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.) 364 DEALLOCATE(lon_ini,lat_ini) 365 366 !--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0 367 WRITE(lunout,*) 368 WRITE(lunout,*)'*** Compute surface geopotential ***' 369 370 !--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS 371 CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque) 372 phis = phis * 9.81 373 phis(iml,:) = phis(1,:) 374 DEALLOCATE(relief_hi,lon_rad,lat_rad) 375 376 END SUBROUTINE start_init_orog0 377 ! 378 !------------------------------------------------------------------------------- 379 380 381 !------------------------------------------------------------------------------- 382 ! 383 SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask) 384 ! 385 !=============================================================================== 386 ! Purpose: Extracted from grid_noro to provide geopotential height for dynamics 387 ! without any call to physics subroutines. 388 !=============================================================================== 389 IMPLICIT NONE 390 !------------------------------------------------------------------------------- 391 ! Arguments: 392 REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp) 393 REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp) 394 REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar) 395 REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar) 396 REAL, INTENT(INOUT):: mask(:,:) !--- MASK (imar+1,jmar) 397 !------------------------------------------------------------------------------- 398 ! Local variables: 399 CHARACTER(LEN=256) :: modname="grid_noro0" 400 REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2) 401 REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2) 402 REAL, ALLOCATABLE :: weight(:,:) ! dim (imar+1,jmar) 403 REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:) ! dim (imar+1,jmar) 404 REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax) 405 REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax) 406 REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax) 407 LOGICAL :: masque_lu 408 INTEGER :: i, ii, imdp, imar, iext 409 INTEGER :: j, jj, jmdp, jmar, nn 410 REAL :: xpi, zlenx, weighx, xincr, zbordnor, zmeanor, zweinor, zbordest 411 REAL :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue 412 !------------------------------------------------------------------------------- 413 imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp") 414 jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp") 415 imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1 416 jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar") 417 IF(imar/=iim) CALL abort_gcm(TRIM(modname),'imar/=iim' ,1) 418 IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1) 419 iext=imdp/10 420 xpi = ACOS(-1.) 421 rad = 6371229. 422 423 !--- ARE WE USING A READ MASK ? 424 masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0 425 WRITE(lunout,*)'Masque lu: ',masque_lu 426 427 !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES: 428 ALLOCATE(xusn(imdp+2*iext)) 429 xusn(1 +iext:imdp +iext)=xd(:) 430 xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi 431 xusn(1+imdp+iext:imdp+2*iext)=xd(1 :iext)+2.*xpi 432 433 ALLOCATE(yusn(jmdp+2)) 434 yusn(1 )=yd(1) +(yd(1) -yd(2)) 435 yusn(2:jmdp+1)=yd(:) 436 yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1)) 437 438 ALLOCATE(zusn(imdp+2*iext,jmdp+2)) 439 zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :) 440 zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :) 441 zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :) 442 zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2) 443 zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2) 444 zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1) 445 zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1) 446 447 !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) 448 ALLOCATE(a(imar+1),b(imar+1)) 449 b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0 450 b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0 451 a(1)=x(1)-(x(2)-x(1))/2.0 452 a(2:imar+1)= b(1:imar) 453 454 ALLOCATE(c(jmar),d(jmar)) 455 d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0 456 d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0 457 c(1)=y(1)-(y(2)-y(1))/2.0 458 c(2:jmar)=d(1:jmar-1) 459 460 !--- INITIALIZATIONS: 461 ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0 462 ALLOCATE(zmea (imar+1,jmar)); zmea (:,:)= 0.0 463 464 !--- SUMMATION OVER GRIDPOINT AREA 465 zleny=xpi/REAL(jmdp)*rad 466 xincr=xpi/REAL(jmdp)/2. 467 ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0. 468 ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0. 469 DO ii = 1, imar+1 470 DO jj = 1, jmar 471 DO j = 2,jmdp+1 472 zlenx =zleny *COS(yusn(j)) 473 zbordnor=(xincr+c(jj)-yusn(j))*rad 474 zbordsud=(xincr-d(jj)+yusn(j))*rad 475 weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny)) 476 IF(weighy/=0) THEN 477 DO i = 2, imdp+2*iext-1 478 zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j)) 479 zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j)) 480 weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx)) 481 IF(weighx/=0)THEN 482 num_tot(ii,jj)=num_tot(ii,jj)+1.0 483 IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0 484 weight(ii,jj)=weight(ii,jj)+weighx*weighy 485 zmea (ii,jj)=zmea (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN 486 END IF 487 END DO 488 END IF 489 END DO 490 END DO 491 END DO 492 493 !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME 494 IF(.NOT.masque_lu) THEN 495 WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:) 496 END IF 497 nn=COUNT(weight(:,1:jmar-1)==0.0) 498 IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn 499 WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:) 500 501 !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS) 502 ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0 503 WHERE(mask>=0.1) mask_tmp = 1. 504 WHERE(weight(:,:)/=0.0) 505 zphi(:,:)=mask_tmp(:,:)*zmea(:,:) 506 zmea(:,:)=mask_tmp(:,:)*zmea(:,:) 507 END WHERE 508 WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea) 509 510 !--- Values at poles 511 zphi(imar+1,:)=zphi(1,:) 512 513 zweinor=SUM(weight(1:imar, 1),DIM=1) 514 zweisud=SUM(weight(1:imar,jmar),DIM=1) 515 zmeanor=SUM(weight(1:imar, 1)*zmea(1:imar, 1),DIM=1) 516 zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1) 517 zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud 518 519 END SUBROUTINE grid_noro0 520 ! 521 !------------------------------------------------------------------------------- 522 523 524 !------------------------------------------------------------------------------- 525 ! 526 SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol) 296 SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol) 527 297 ! 528 298 !------------------------------------------------------------------------------- … … 534 304 REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) 535 305 REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) 536 LOGICAL, INTENT(IN) :: ibar537 306 REAL, INTENT(IN) :: zs (:,:) ! dim (iml,jml) 538 307 REAL, INTENT(OUT) :: psol(:,:) ! dim (iml,jml) … … 606 375 ALLOCATE(field(iml,jml)) 607 376 CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana) 608 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)609 CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,&610 lon_in, lat_in, lon_in2, lat_in2, field)377 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.) 378 CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana, & 379 lon_in, lat_in, lon_in2, lat_in2, field) 611 380 ELSE IF(SIZE(field)/=SIZE(z)) THEN 612 381 msg='The '//TRIM(msg)//' field we have does not have the right size' … … 625 394 !------------------------------------------------------------------------------- 626 395 ! 627 SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d ,ibar)396 SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d) 628 397 ! 629 398 !------------------------------------------------------------------------------- … … 639 408 REAL, INTENT(IN) :: pls_in(:,:,:) ! dim (iml,jml,lml) 640 409 REAL, INTENT(OUT) :: var3d (:,:,:) ! dim (iml,jml,lml) 641 LOGICAL, INTENT(IN) :: ibar642 410 !------------------------------------------------------------------------------- 643 411 ! Local variables: … … 667 435 ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn)) 668 436 CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, & 669 lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)437 lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.) 670 438 DEALLOCATE(lon_ini, lat_ini) 671 439 … … 673 441 ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) 674 442 DO il = 1,llm_dyn 675 CALL interp_startvar(var,i bar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),&443 CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il), & 676 444 lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il)) 677 445 END DO … … 708 476 !------------------------------------------------------------------------------- 709 477 ! 710 SUBROUTINE interp_startvar(nam,ib ar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)478 SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo) 711 479 ! 712 480 !------------------------------------------------------------------------------- 713 481 USE inter_barxy_m, ONLY: inter_barxy 714 USE grid_atob_m, ONLY: grille_m715 482 IMPLICIT NONE 716 483 !------------------------------------------------------------------------------- 717 484 ! Arguments: 718 485 CHARACTER(LEN=*), INTENT(IN) :: nam 719 LOGICAL, INTENT(IN) :: ib ar, ibeg486 LOGICAL, INTENT(IN) :: ibeg 720 487 REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) 721 488 REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) … … 735 502 j2=SIZE(lat2) 736 503 ALLOCATE(vtmp(i1-1,j1)) 737 IF(ibar) THEN 738 IF(ibeg.AND.prt_level>1) THEN 739 WRITE(lunout,*)"---------------------------------------------------------" 740 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 741 WRITE(lunout,*)"---------------------------------------------------------" 742 END IF 743 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 744 ELSE 745 CALL grille_m (lon, lat, vari, lon1, lat1, vtmp) 746 END IF 504 IF(ibeg.AND.prt_level>1) THEN 505 WRITE(lunout,*)"---------------------------------------------------------" 506 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 507 WRITE(lunout,*)"---------------------------------------------------------" 508 END IF 509 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 747 510 CALL gr_int_dyn(vtmp, varo, i1-1, j1) 748 511 -
LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90
r2320 r2336 36 36 USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo 37 37 USE assert_eq_m, ONLY: assert_eq 38 USE dimphy 38 USE dimphy, ONLY: klon 39 USE conf_dat_m, ONLY: conf_dat2d 39 40 USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, & 40 41 rlon, solsw, radsol, t_ancien, wake_deltat, wake_s, rain_fall, qsol, z0h, & … … 51 52 include "paramet.h" 52 53 include "comgeom2.h" 53 include "comvert.h"54 54 include "comconst.h" 55 55 include "dimsoil.h" 56 56 include "temps.h" 57 include "comdissnew.h"58 include "serre.h"59 57 include "clesphys.h" 60 58 REAL, SAVE :: deg2rad 61 59 REAL, SAVE, ALLOCATABLE :: tsol(:) 62 ! REAL, SAVE, ALLOCATABLE :: rugo(:,:) ! ??? COMPUTED BUT NOT USED ???63 60 INTEGER, SAVE :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys 64 61 REAL, ALLOCATABLE, SAVE :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:) 65 CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc" 66 CHARACTER(LEN=256), PARAMETER :: title="RELIEF" 62 CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF" 63 CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc", psrfvar="SP" 64 CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW", tsrfvar="ST" 67 65 68 66 … … 72 70 !------------------------------------------------------------------------------- 73 71 ! 74 SUBROUTINE etat0phys_netcdf( ib,masque, phis)72 SUBROUTINE etat0phys_netcdf(masque, phis) 75 73 ! 76 74 !------------------------------------------------------------------------------- 77 75 ! Purpose: Creates initial states 78 76 !------------------------------------------------------------------------------- 79 ! Note: This routine is designed to work for Earth 77 ! Notes: 1) This routine is designed to work for Earth 78 ! 2) If masque(:,:)/=-99999., masque and phis are already known. 79 ! Otherwise: compute it. 80 80 !------------------------------------------------------------------------------- 81 81 USE control_mod … … 84 84 USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz 85 85 USE indice_sol_mod 86 USE conf_phys_m, ONLY: conf_phys 87 USE exner_hyb_m, ONLY: exner_hyb 88 USE exner_milieu_m, ONLY: exner_milieu 89 USE test_disvert_m, ONLY: test_disvert 90 USE grid_atob_m, ONLY: grille_m 86 USE conf_phys_m, ONLY: conf_phys 87 USE init_ssrf_m, ONLY: start_init_subsurf 91 88 IMPLICIT NONE 92 89 !------------------------------------------------------------------------------- 93 90 ! Arguments: 94 LOGICAL, INTENT(IN) :: ib !--- Barycentric interpolation95 91 REAL, INTENT(INOUT) :: masque(:,:) !--- Land mask dim(iip1,jjp1) 96 92 REAL, INTENT(INOUT) :: phis (:,:) !--- Ground geopotential dim(iip1,jjp1) … … 99 95 CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt 100 96 INTEGER :: i, j, l, ji, iml, jml 101 REAL :: phystep 97 LOGICAL :: read_mask 98 REAL :: phystep, dummy 102 99 REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp 103 100 REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0, fder 104 101 REAL, DIMENSION(klon,nbsrf) :: qsolsrf, snsrf 105 102 REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil 106 107 !--- Local variables for sea-ice reading:108 LOGICAL :: read_mask109 INTEGER :: iml_lic, jml_lic, isst(klon-2)110 INTEGER :: fid, llm_tmp, ttm_tmp, itaul(1)111 REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:)112 REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:)113 REAL :: date, lev(1), dummy114 REAL :: flic_tmp(SIZE(masque,1),SIZE(masque,2))115 103 116 104 !--- Arguments for conf_phys … … 132 120 jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml") 133 121 134 ! Grid construction and miscellanous initializations.122 ! Physics configuration 135 123 !******************************************************************************* 136 124 CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & … … 145 133 read_climoz, & 146 134 alp_offset) 147 148 135 CALL phys_state_var_init(read_climoz) 149 136 150 co2_ppm0 = co2_ppm !--- Initial atmospheric CO2 conc. from .def file 151 152 rlat(1) = pi/2. 153 DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j); END DO 154 rlat(klon) = - pi/2. 155 rlat(:)=rlat(:)*(180.0/pi) 156 157 rlon(1) = 0.0 158 DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO 159 rlon(klon) = 0.0 160 rlon(:)=rlon(:)*(180.0/pi) 137 !--- Initial atmospheric CO2 conc. from .def file 138 co2_ppm0 = co2_ppm 161 139 162 140 ! Compute ground geopotential, sub-cells quantities and possibly the mask. … … 174 152 WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1. 175 153 END IF 176 CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid154 CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid 177 155 178 156 ! Compute tsol and qsol on physical grid, knowing phis on 2D grid. 179 157 !******************************************************************************* 180 CALL start_init_phys(rlon v, rlatu, rlonu, rlatv, ib, phis)158 CALL start_init_phys(rlonu, rlatv, phis) 181 159 182 160 ! Some initializations. … … 188 166 CALL regr_lat_time_climoz(read_climoz) 189 167 190 ! Sub-surfaces initialization 191 !******************************************************************************* 192 !--- Read and interpolate on model T-grid soil fraction and soil ice fraction. 193 CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid) 194 ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic)) 195 ALLOCATE(dlat_lic(jml_lic), dlon_lic(iml_lic)) 196 ALLOCATE( fraclic(iml_lic,jml_lic)) 197 CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp, & 198 & lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) 199 CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic) 200 CALL flinclo(fid) 201 WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic 202 IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees 203 IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180. 204 dlon_lic(:)=lon_lic(:,1) 205 dlat_lic(:)=lat_lic(1,:) 206 CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) ) 207 flic_tmp(iml,:)=flic_tmp(1,:) 208 209 !--- To the physical grid 210 pctsrf(:,:) = 0. 211 CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic)) 212 213 !--- Adequation with soil/sea mask 214 WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 215 WHERE(zmasq(:)<EPSFRA) pctsrf(:,is_lic)=0. 216 pctsrf(:,is_ter)=zmasq(:) 217 DO ji=1,klon 218 IF(zmasq(ji)>EPSFRA) THEN 219 IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN 220 pctsrf(ji,is_lic)=zmasq(ji) 221 pctsrf(ji,is_ter)=0. 222 ELSE 223 pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic) 224 IF(pctsrf(ji,is_ter)<EPSFRA) THEN 225 pctsrf(ji,is_ter)=0. 226 pctsrf(ji,is_lic)=zmasq(ji) 227 END IF 228 END IF 229 END IF 230 END DO 231 232 !--- Sub-surface ocean and sea ice (sea ice set to zero for start). 233 pctsrf(:,is_oce)=(1.-zmasq(:)) 234 WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0. 235 IF(read_mask) pctsrf(:,is_oce)=1-zmasq(:) 236 isst=0 237 WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1 238 239 !--- It is checked that the sub-surfaces sum is equal to 1. 240 ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA) 241 IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points' 168 ! Sub-surfaces initialization. 169 !******************************************************************************* 170 CALL start_init_subsurf(read_mask) 242 171 243 172 ! Write physical initial state … … 321 250 ! described in LOTT & MILLER (1997) and LOTT(1999). 322 251 !=============================================================================== 323 USE conf_dat_m, ONLY: conf_dat2d324 ! USE grid_atob_m, ONLY: rugsoro325 252 USE grid_noro_m, ONLY: grid_noro 326 253 IMPLICIT NONE … … 331 258 !------------------------------------------------------------------------------- 332 259 ! Local variables: 333 CHARACTER(LEN=256) :: modname="start_init_orog" 334 CHARACTER(LEN=256) :: title="RELIEF" 260 CHARACTER(LEN=256) :: modname 335 261 INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1) 336 262 REAL :: lev(1), date, dt … … 340 266 REAL, ALLOCATABLE :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:) 341 267 !------------------------------------------------------------------------------- 268 modname="start_init_orog" 342 269 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml") 343 270 jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml") … … 350 277 lev, ttm_tmp, itau, date, dt, fid) 351 278 ALLOCATE(relief_hi(iml_rel,jml_rel)) 352 CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)279 CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi) 353 280 CALL flinclo(fid) 354 281 … … 360 287 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 361 288 ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) 362 CALL conf_dat2d( title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)289 CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.) 363 290 DEALLOCATE(lon_ini,lat_ini) 364 291 … … 378 305 phis = phis * 9.81 379 306 phis(iml,:) = phis(1,:) 380 381 !--- COMPUTE SURFACE ROUGHNESS382 ! WRITE(lunout,*)383 ! WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***'384 ! ALLOCATE(tmp_var(iml-1,jml))385 ! CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var)386 ! ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)387 ! DEALLOCATE(tmp_var)388 307 DEALLOCATE(relief_hi,lon_rad,lat_rad) 389 308 … … 405 324 !------------------------------------------------------------------------------- 406 325 ! 407 SUBROUTINE start_init_phys(lon_in,lat_in, lon_in2,lat_in2,ibar,phis)326 SUBROUTINE start_init_phys(lon_in,lat_in,phis) 408 327 ! 409 328 !=============================================================================== … … 413 332 !------------------------------------------------------------------------------- 414 333 ! Arguments: 415 REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) 416 REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) 417 LOGICAL, INTENT(IN) :: ibar 334 REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml2) 418 335 REAL, INTENT(IN) :: phis(:,:) ! dim (iml,jml) 419 336 !------------------------------------------------------------------------------- 420 337 ! Local variables: 421 CHARACTER(LEN=256) :: modname ="start_init_phys", physfname="ECPHY.nc"338 CHARACTER(LEN=256) :: modname 422 339 REAL :: date, dt 423 340 INTEGER :: iml, jml, jml2, itau(1) … … 426 343 REAL, ALLOCATABLE :: ts(:,:), qs(:,:) 427 344 !------------------------------------------------------------------------------- 428 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml")429 jml=assert_eq(SIZE(lat_in),SIZE(phis,2), TRIM(modname)//" jml")430 jml 2=SIZE(lat_in2)345 modname="start_init_phys" 346 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml") 347 jml=SIZE(phis,2); jml2=SIZE(lat_in) 431 348 432 349 WRITE(lunout,*)'Opening the surface analysis' 433 CALL flininfo(phy sfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)434 WRITE(lunout,*) 'Values read: ', 350 CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys) 351 WRITE(lunout,*) 'Values read: ', iml_phys, jml_phys, llm_phys, ttm_phys 435 352 436 353 ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys)) 437 354 ALLOCATE(levphys_ini(llm_phys)) 438 CALL flinopen(phy sfname, .FALSE., iml_phys, jml_phys, llm_phys, &355 CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys, & 439 356 lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys) 440 357 … … 445 362 446 363 ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys)) 447 CALL get_var_phys( 'ST',ts) !--- SURFACE TEMPERATURE448 CALL get_var_phys( 'CDSW',qs) !--- SOIL MOISTURE364 CALL get_var_phys(tsrfvar,ts) !--- SURFACE TEMPERATURE 365 CALL get_var_phys(qsolvar,qs) !--- SOIL MOISTURE 449 366 CALL flinclo(fid_phys) 450 367 DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) … … 463 380 ! 464 381 !------------------------------------------------------------------------------- 465 USE conf_dat_m, ONLY: conf_dat2d466 382 IMPLICIT NONE 467 383 !------------------------------------------------------------------------------- … … 474 390 !------------------------------------------------------------------------------- 475 391 SELECT CASE(title) 476 CASE( 'SP');tllm=0477 CASE( 'ST','CDSW'); tllm=llm_phys392 CASE(psrfvar); tllm=0 393 CASE(tsrfvar,qsolvar); tllm=llm_phys 478 394 END SELECT 479 395 IF(ALLOCATED(field)) RETURN 480 396 ALLOCATE(field(iml,jml)); field(:,:)=0. 481 397 CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana) 482 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)483 CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,&484 lon_in, lat_in, lon_in2, lat_in2,field)398 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.) 399 CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana, & 400 lon_in, lat_in, field) 485 401 486 402 END SUBROUTINE get_var_phys … … 495 411 !------------------------------------------------------------------------------- 496 412 ! 497 SUBROUTINE interp_startvar(nam,ib ar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)413 SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo) 498 414 ! 499 415 !------------------------------------------------------------------------------- 500 416 USE inter_barxy_m, ONLY: inter_barxy 501 USE grid_atob_m, ONLY: grille_m502 417 IMPLICIT NONE 503 418 !------------------------------------------------------------------------------- 504 419 ! Arguments: 505 420 CHARACTER(LEN=*), INTENT(IN) :: nam 506 LOGICAL, INTENT(IN) :: ib ar, ibeg421 LOGICAL, INTENT(IN) :: ibeg 507 422 REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) 508 423 REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) 509 REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1)510 424 REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) 511 425 REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1) 512 426 !------------------------------------------------------------------------------- 513 427 ! Local variables: 514 CHARACTER(LEN=256) :: modname ="interp_startvar"428 CHARACTER(LEN=256) :: modname 515 429 INTEGER :: ii, jj, i1, j1, j2 516 430 REAL, ALLOCATABLE :: vtmp(:,:) 517 431 !------------------------------------------------------------------------------- 518 ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")519 jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")520 i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")521 j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1")522 j 2=SIZE(lat2)432 modname="interp_startvar" 433 ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii") 434 jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj") 435 i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1") 436 j1=SIZE(varo,2); j2=SIZE(lat2) 523 437 ALLOCATE(vtmp(i1-1,j1)) 524 IF(ibar) THEN 525 IF(ibeg.AND.prt_level>1) THEN 526 WRITE(lunout,*)"--------------------------------------------------------" 527 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 528 WRITE(lunout,*)"--------------------------------------------------------" 529 END IF 530 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 531 ELSE 532 CALL grille_m (lon, lat, vari, lon1, lat1, vtmp) 438 IF(ibeg.AND.prt_level>1) THEN 439 WRITE(lunout,*)"--------------------------------------------------------" 440 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 441 WRITE(lunout,*)"--------------------------------------------------------" 533 442 END IF 443 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 534 444 CALL gr_int_dyn(vtmp, varo, i1-1, j1) 535 445 -
LMDZ5/trunk/libf/phylmd/limit_netcdf.F90
r2315 r2336 1 SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque) 2 ! 3 ! -------------------------------------------------------------------------------1 MODULE limit 2 ! 3 !******************************************************************************* 4 4 ! Author : L. Fairhead, 27/01/94 5 5 !------------------------------------------------------------------------------- … … 16 16 ! * 12/2009: D. Cugnet (f77->f90, calendars, files from coupled runs) 17 17 !------------------------------------------------------------------------------- 18 #ifndef CPP_1D 18 19 USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, & 20 ioconf_calendar, ioget_calendar, lock_calendar, ioget_mon_len, ioget_year_len 21 USE assert_eq_m, ONLY: assert_eq 22 USE conf_dat_m, ONLY: conf_dat2d, conf_dat3d 23 USE dimphy, ONLY: klon, zmasq 24 USE phys_state_var_mod, ONLY: pctsrf, rlon, rlat 19 25 USE control_mod 20 USE indice_sol_mod 21 USE dimphy 22 USE ioipsl, ONLY : ioget_year_len 23 USE phys_state_var_mod, ONLY : pctsrf, rlon, rlat 24 USE netcdf, ONLY : NF90_OPEN, NF90_CREATE, NF90_CLOSE, & 25 NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT, & 26 NF90_NOERR, NF90_NOWRITE, NF90_DOUBLE, NF90_GLOBAL, & 27 NF90_CLOBBER, NF90_ENDDEF, NF90_UNLIMITED, NF90_FLOAT 28 USE inter_barxy_m, only: inter_barxy 29 USE netcdf95, ONLY: nf95_def_var, nf95_put_att, nf95_put_var 30 USE grid_atob_m, ONLY: grille_m, rugosite, sea_ice 31 USE print_control_mod, ONLY: prt_level,lunout 32 IMPLICIT NONE 33 !------------------------------------------------------------------------------- 34 ! Arguments: 35 include "dimensions.h" 36 include "paramet.h" 37 LOGICAL, INTENT(IN) :: interbar ! barycentric interpolation 38 LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag 39 LOGICAL, INTENT(IN) :: oldice ! old way ice computation 40 REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask 41 !------------------------------------------------------------------------------- 42 ! Local variables: 43 include "logic.h" 44 include "comgeom2.h" 45 include "comconst.h" 46 47 !--- INPUT NETCDF FILES NAMES -------------------------------------------------- 48 CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam 26 USE init_ssrf_m, ONLY: start_init_subsurf 27 49 28 CHARACTER(LEN=20), PARAMETER :: & 50 29 fsst(4)=['amipbc_sst_1x1.nc ','cpl_atm_sst.nc ','histmth_sst.nc '& … … 56 35 vsst(4)=['tosbcs ','SISUTESW ','tsol_oce ','sstk '], & 57 36 vsic(4)=['sicbcs ','SIICECOV ','pourc_sic ','ci '] 58 CHARACTER(LEN=20), PARAMETER :: frugo='Rugos.nc ', & 59 falbe='Albedo.nc ' 60 CHARACTER(LEN=10), PARAMETER :: vrug='RUGOS ', valb='ALBEDO ' 37 CHARACTER(LEN=10), PARAMETER :: & 38 frugo='Rugos.nc ', falbe='Albedo.nc ', frelf='Relief.nc ', & 39 vrug='RUGOS ', valb='ALBEDO ', vrel='RELIEF ' 40 41 CONTAINS 42 43 !------------------------------------------------------------------------------- 44 ! 45 SUBROUTINE limit_netcdf(masque, phis, extrap) 46 ! 47 !------------------------------------------------------------------------------- 48 ! Author : L. Fairhead, 27/01/94 49 !------------------------------------------------------------------------------- 50 ! Purpose: Boundary conditions files building for new model using climatologies. 51 ! Both grids have to be regular. 52 !------------------------------------------------------------------------------- 53 ! Note: This routine is designed to work for Earth 54 !------------------------------------------------------------------------------- 55 ! Modification history: 56 ! * 23/03/1994: Z. X. Li 57 ! * 09/1999: L. Fairhead (netcdf reading in LMDZ.3.3) 58 ! * 07/2001: P. Le Van 59 ! * 11/2009: L. Guez (ozone day & night climatos, see etat0_netcdf.F90) 60 ! * 12/2009: D. Cugnet (f77->f90, calendars, files from coupled runs) 61 !------------------------------------------------------------------------------- 62 #ifndef CPP_1D 63 USE indice_sol_mod 64 USE netcdf, ONLY: NF90_OPEN, NF90_CREATE, NF90_CLOSE, & 65 NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT, & 66 NF90_NOERR, NF90_NOWRITE, NF90_DOUBLE, NF90_GLOBAL, & 67 NF90_CLOBBER, NF90_ENDDEF, NF90_UNLIMITED, NF90_FLOAT 68 USE inter_barxy_m, ONLY: inter_barxy 69 USE netcdf95, ONLY: nf95_def_var, nf95_put_att, nf95_put_var 70 IMPLICIT NONE 71 !------------------------------------------------------------------------------- 72 ! Arguments: 73 include "iniprint.h" 74 include "dimensions.h" 75 include "paramet.h" 76 REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask 77 REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: phis ! ground geopotential 78 LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag 79 !------------------------------------------------------------------------------- 80 ! Local variables: 81 include "logic.h" 82 include "comgeom2.h" 83 include "comconst.h" 84 85 !--- INPUT NETCDF FILES NAMES -------------------------------------------------- 86 CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam 61 87 CHARACTER(LEN=10) :: varname 88 62 89 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------ 63 90 REAL :: fi_ice(klon), verif(klon) … … 82 109 CALL inigeom 83 110 111 !--- MASK, GROUND GEOPOT. & SUBSURFACES COMPUTATION (IN CASE ok_etat0==.FALSE.) 112 IF(ALL(masque==-99999.)) THEN 113 CALL start_init_orog0(rlonv,rlatu,phis,masque) 114 CALL gr_dyn_fi(1,iip1,jjp1,klon,masque,zmasq) !--- To physical grid 115 ALLOCATE(rlon(klon),rlat(klon),pctsrf(klon,nbsrf)) 116 CALL start_init_subsurf(.FALSE.) 117 END IF 118 84 119 !--- Beware: anneeref (from gcm.def) is used to determine output time sampling 85 120 ndays=ioget_year_len(anneeref) … … 87 122 !--- RUGOSITY TREATMENT -------------------------------------------------------- 88 123 CALL msg(1,'Traitement de la rugosite') 89 CALL get_2Dfield(frugo,vrug,'RUG', interbar,ndays,phy_rug,mask=masque(1:iim,:))124 CALL get_2Dfield(frugo,vrug,'RUG',ndays,phy_rug,mask=masque(1:iim,:)) 90 125 91 126 !--- OCEAN TREATMENT ----------------------------------------------------------- … … 108 143 CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile)) 109 144 110 CALL get_2Dfield(icefile,varname, 'SIC', interbar,ndays,phy_ice,flag=oldice)145 CALL get_2Dfield(icefile,varname, 'SIC',ndays,phy_ice) 111 146 112 147 ALLOCATE(pctsrf_t(klon,nbsrf,ndays)) … … 167 202 CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile)) 168 203 169 CALL get_2Dfield(sstfile,varname,'SST', interbar,ndays,phy_sst,flag=extrap)204 CALL get_2Dfield(sstfile,varname,'SST',ndays,phy_sst,flag=extrap) 170 205 171 206 !--- ALBEDO TREATMENT ---------------------------------------------------------- 172 207 CALL msg(1,'Traitement de l albedo') 173 CALL get_2Dfield(falbe,valb,'ALB', interbar,ndays,phy_alb)208 CALL get_2Dfield(falbe,valb,'ALB',ndays,phy_alb) 174 209 175 210 !--- REFERENCE GROUND HEAT FLUX TREATMENT -------------------------------------- … … 251 286 !------------------------------------------------------------------------------- 252 287 ! 253 SUBROUTINE get_2Dfield(fnam, varname, mode, ibar,ndays, champo, flag, mask)288 SUBROUTINE get_2Dfield(fnam, varname, mode, ndays, champo, flag, mask) 254 289 ! 255 290 !----------------------------------------------------------------------------- … … 263 298 NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, & 264 299 NF90_GET_ATT 265 USE dimphy, ONLY : klon266 USE phys_state_var_mod, ONLY : pctsrf267 USE conf_dat_m, ONLY: conf_dat2d268 USE control_mod269 300 USE pchsp_95_m, only: pchsp_95 270 301 USE pchfe_95_m, only: pchfe_95 … … 281 312 CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name 282 313 CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB 283 LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels284 314 INTEGER, INTENT(IN) :: ndays ! current year number of days 285 315 REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) … … 312 342 CHARACTER(LEN=25) :: title ! for messages 313 343 LOGICAL :: extrp ! flag for extrapolation 314 LOGICAL :: oldice ! flag for old way ice computation315 344 REAL :: chmin, chmax 316 345 INTEGER ierr … … 328 357 CASE('ALB'); title='Albedo' 329 358 END SELECT 330 331 332 extrp=.FALSE. 333 oldice=.FALSE. 334 IF ( PRESENT(flag) ) THEN 335 IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 336 IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 337 END IF 359 extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag 338 360 339 361 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- … … 400 422 DO l=1, lmdep 401 423 CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam) 402 CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, ibar)424 CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.) 403 425 IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work) 404 405 IF(ibar .AND. .NOT.oldice) THEN 406 IF(l==1) THEN 426 IF(l==1) THEN 407 427 CALL msg(5,"----------------------------------------------------------") 408 428 CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$") 409 429 CALL msg(5,"----------------------------------------------------------") 410 END IF 411 IF(mode=='RUG') champ=LOG(champ) 412 CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint) 413 IF(mode=='RUG') THEN 414 champint=EXP(champint) 415 WHERE(NINT(mask)/=1) champint=0.001 416 END IF 417 ELSE 418 SELECT CASE(mode) 419 CASE('RUG'); CALL rugosite(dlon,dlat,champ,rlonv,rlatu,champint,mask) 420 CASE('SIC'); CALL sea_ice (dlon,dlat,champ,rlonv,rlatu,champint) 421 CASE('SST','ALB'); CALL grille_m(dlon,dlat,champ,rlonv,rlatu,champint) 422 END SELECT 430 END IF 431 IF(mode=='RUG') champ=LOG(champ) 432 CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint) 433 IF(mode=='RUG') THEN 434 champint=EXP(champint) 435 WHERE(NINT(mask)/=1) champint=0.001 423 436 END IF 424 437 champtime(:, :, l)=champint … … 501 514 !------------------------------------------------------------------------------- 502 515 ! 516 SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque) 517 ! 518 !------------------------------------------------------------------------------- 519 IMPLICIT NONE 520 !=============================================================================== 521 ! Purpose: Compute "phis" just like it would be in start_init_orog. 522 !=============================================================================== 523 ! Arguments: 524 REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) 525 REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml) 526 !------------------------------------------------------------------------------- 527 ! Local variables: 528 CHARACTER(LEN=256) :: modname="start_init_orog0" 529 INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1) 530 REAL :: lev(1), date, dt, deg2rad 531 REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:) 532 REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:) 533 !------------------------------------------------------------------------------- 534 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml") 535 jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml") 536 IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1) 537 IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1) 538 pi=2.0*ASIN(1.0); deg2rad=pi/180.0 539 IF(ANY(phis/=-99999.)) RETURN !--- phis ALREADY KNOWN 540 541 !--- HIGH RESOLUTION OROGRAPHY 542 CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) 543 544 ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) 545 CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel, & 546 lev, ttm_tmp, itau, date, dt, fid) 547 ALLOCATE(relief_hi(iml_rel,jml_rel)) 548 CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) 549 CALL flinclo(fid) 550 551 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 552 ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) 553 lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad 554 lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad 555 556 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 557 ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) 558 CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.) 559 DEALLOCATE(lon_ini,lat_ini) 560 561 !--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0 562 WRITE(lunout,*) 563 WRITE(lunout,*)'*** Compute surface geopotential ***' 564 565 !--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS 566 CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque) 567 phis = phis * 9.81 568 phis(iml,:) = phis(1,:) 569 DEALLOCATE(relief_hi,lon_rad,lat_rad) 570 571 END SUBROUTINE start_init_orog0 572 ! 573 !------------------------------------------------------------------------------- 574 575 576 !------------------------------------------------------------------------------- 577 ! 578 SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask) 579 ! 580 !=============================================================================== 581 ! Purpose: Extracted from grid_noro to provide geopotential height for dynamics 582 ! without any call to physics subroutines. 583 !=============================================================================== 584 IMPLICIT NONE 585 !------------------------------------------------------------------------------- 586 ! Arguments: 587 REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp) 588 REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp) 589 REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar) 590 REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar) 591 REAL, INTENT(INOUT):: mask(:,:) !--- MASK (imar+1,jmar) 592 !------------------------------------------------------------------------------- 593 ! Local variables: 594 CHARACTER(LEN=256) :: modname="grid_noro0" 595 REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2) 596 REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2) 597 REAL, ALLOCATABLE :: weight(:,:) ! dim (imar+1,jmar) 598 REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:) ! dim (imar+1,jmar) 599 REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax) 600 REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax) 601 REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax) 602 LOGICAL :: masque_lu 603 INTEGER :: i, ii, imdp, imar, iext 604 INTEGER :: j, jj, jmdp, jmar, nn 605 REAL :: xpi, zlenx, weighx, xincr, zbordnor, zmeanor, zweinor, zbordest 606 REAL :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue 607 !------------------------------------------------------------------------------- 608 imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp") 609 jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp") 610 imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1 611 jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar") 612 IF(imar/=iim) CALL abort_gcm(TRIM(modname),'imar/=iim' ,1) 613 IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1) 614 iext=imdp/10 615 xpi = ACOS(-1.) 616 rad = 6371229. 617 618 !--- ARE WE USING A READ MASK ? 619 masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0 620 WRITE(lunout,*)'Masque lu: ',masque_lu 621 622 !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES: 623 ALLOCATE(xusn(imdp+2*iext)) 624 xusn(1 +iext:imdp +iext)=xd(:) 625 xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi 626 xusn(1+imdp+iext:imdp+2*iext)=xd(1 :iext)+2.*xpi 627 628 ALLOCATE(yusn(jmdp+2)) 629 yusn(1 )=yd(1) +(yd(1) -yd(2)) 630 yusn(2:jmdp+1)=yd(:) 631 yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1)) 632 633 ALLOCATE(zusn(imdp+2*iext,jmdp+2)) 634 zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :) 635 zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :) 636 zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :) 637 zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2) 638 zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2) 639 zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1) 640 zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1) 641 642 !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) 643 ALLOCATE(a(imar+1),b(imar+1)) 644 b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0 645 b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0 646 a(1)=x(1)-(x(2)-x(1))/2.0 647 a(2:imar+1)= b(1:imar) 648 649 ALLOCATE(c(jmar),d(jmar)) 650 d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0 651 d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0 652 c(1)=y(1)-(y(2)-y(1))/2.0 653 c(2:jmar)=d(1:jmar-1) 654 655 !--- INITIALIZATIONS: 656 ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0 657 ALLOCATE(zmea (imar+1,jmar)); zmea (:,:)= 0.0 658 659 !--- SUMMATION OVER GRIDPOINT AREA 660 zleny=xpi/REAL(jmdp)*rad 661 xincr=xpi/REAL(jmdp)/2. 662 ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0. 663 ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0. 664 DO ii = 1, imar+1 665 DO jj = 1, jmar 666 DO j = 2,jmdp+1 667 zlenx =zleny *COS(yusn(j)) 668 zbordnor=(xincr+c(jj)-yusn(j))*rad 669 zbordsud=(xincr-d(jj)+yusn(j))*rad 670 weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny)) 671 IF(weighy/=0) THEN 672 DO i = 2, imdp+2*iext-1 673 zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j)) 674 zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j)) 675 weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx)) 676 IF(weighx/=0)THEN 677 num_tot(ii,jj)=num_tot(ii,jj)+1.0 678 IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0 679 weight(ii,jj)=weight(ii,jj)+weighx*weighy 680 zmea (ii,jj)=zmea (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN 681 END IF 682 END DO 683 END IF 684 END DO 685 END DO 686 END DO 687 688 !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME 689 IF(.NOT.masque_lu) THEN 690 WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:) 691 END IF 692 nn=COUNT(weight(:,1:jmar-1)==0.0) 693 IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn 694 WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:) 695 696 !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS) 697 ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0 698 WHERE(mask>=0.1) mask_tmp = 1. 699 WHERE(weight(:,:)/=0.0) 700 zphi(:,:)=mask_tmp(:,:)*zmea(:,:) 701 zmea(:,:)=mask_tmp(:,:)*zmea(:,:) 702 END WHERE 703 WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea) 704 705 !--- Values at poles 706 zphi(imar+1,:)=zphi(1,:) 707 708 zweinor=SUM(weight(1:imar, 1),DIM=1) 709 zweisud=SUM(weight(1:imar,jmar),DIM=1) 710 zmeanor=SUM(weight(1:imar, 1)*zmea(1:imar, 1),DIM=1) 711 zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1) 712 zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud 713 714 END SUBROUTINE grid_noro0 715 ! 716 !------------------------------------------------------------------------------- 717 718 719 !------------------------------------------------------------------------------- 720 ! 503 721 FUNCTION year_len(y,cal_in) 504 722 ! 505 723 !------------------------------------------------------------------------------- 506 USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len507 724 IMPLICIT NONE 508 725 !------------------------------------------------------------------------------- … … 537 754 ! 538 755 !------------------------------------------------------------------------------- 539 USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len540 756 IMPLICIT NONE 541 757 !------------------------------------------------------------------------------- … … 624 840 !------------------------------------------------------------------------------- 625 841 USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR 626 USE print_control_mod, ONLY: lunout627 842 IMPLICIT NONE 628 843 !------------------------------------------------------------------------------- … … 643 858 ! of #ifndef CPP_1D 644 859 END SUBROUTINE limit_netcdf 860 861 END MODULE limit 862 ! 863 !******************************************************************************* -
LMDZ5/trunk/makelmdz_fcm
r2326 r2336 184 184 done 185 185 186 if [[ $code = ce0l && $paramem = mem ]]187 then188 echo "There is no parallel version of ce0l at the moment."189 echo "Please compile the sequential version of the code to produce the" \190 "executable ce0l."191 exit 1192 fi193 194 186 ############################################################### 195 187 # path to fcm
Note: See TracChangeset
for help on using the changeset viewer.
