Changeset 4206 for dynamico_lmdz/simple_physics/phyparam
- Timestamp:
- Dec 22, 2019, 1:10:51 AM (5 years ago)
- Location:
- dynamico_lmdz/simple_physics/phyparam/physics
- Files:
-
- 3 edited
-
convection.F90 (modified) (3 diffs)
-
logging.F90 (modified) (5 diffs)
-
use_logging.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
dynamico_lmdz/simple_physics/phyparam/physics/convection.F90
r4205 r4206 1 1 MODULE convection 2 3 #include "use_logging.h" 4 2 5 IMPLICIT NONE 3 6 4 7 CONTAINS 8 9 PURE SUBROUTINE sigma_levels(ngrid, nlay, i, pplev, ppopsk, sig, dsig, sdsig) 10 INTEGER, INTENT(IN) :: ngrid, nlay, i 11 REAL, INTENT(IN) :: pplev(ngrid,nlay+1), ppopsk(ngrid,nlay) 12 REAL, INTENT(OUT) :: sig(nlay+1), dsig(nlay), sdsig(nlay) 13 ! Calcul des niveaux sigma sur cette colonne 14 INTEGER :: l 15 DO l=1,nlay+1 16 sig(l)=pplev(i,l)/pplev(i,1) 17 ENDDO 18 DO l=1,nlay 19 dsig(l)=sig(l)-sig(l+1) 20 sdsig(l)=ppopsk(i,l)*dsig(l) 21 ENDDO 22 END SUBROUTINE sigma_levels 23 24 SUBROUTINE adjust_onecolumn(ngrid, nlay, i, sig, dsig, sdsig, zu2, zv2, zh2) 25 INTEGER, INTENT(IN) :: ngrid, nlay, i 26 REAL, INTENT(OUT) :: sig(nlay+1), sdsig(nlay), dsig(nlay) 27 REAL, INTENT(INOUT) :: zu2(ngrid,nlay), zv2(ngrid,nlay), zh2(ngrid,nlay) 28 29 INTEGER :: l,l1,l2,jj 30 LOGICAL :: extend 31 REAL :: zhm,zsm,zum,zvm,zalpha 32 #include "dimensions.h" 33 34 l2=1 35 DO WHILE( l2<nlay ) 36 l2 = l2+1 37 ! starting from l1=l2, we shall extend l1..l2 downwards and upwards 38 ! until zhm = potential temperature mass-weighted over l1..l2 be 39 ! higher than pot. temp. at level l1-1 40 ! and lower than pot. temp. at level l2+1 41 ! for this we incrementally compute zsm = mass of layers l1 .. l2 and zhm 42 zsm = sdsig(l2) 43 zhm = zh2(i, l2) 44 l1 = l2 45 extend = .TRUE. 46 DO WHILE(extend) 47 ! extend downwards if l1>1 and zhm is lower than layer l1-1 48 extend = .FALSE. 49 IF (l1 > 1) THEN 50 IF (zhm < zh2(i, l1-1)) extend = .TRUE. 51 END IF 52 IF(extend) THEN 53 l1 = l1-1 54 zsm = zsm + sdsig(l1) 55 zhm = zhm + sdsig(l1) * (zh2(i, l1) - zhm) / zsm 56 CYCLE 57 END IF 58 ! extend upwards if l2<nlay and zhm is higher than layer l2+1 59 extend=.FALSE. 60 IF (l2 < nlay) THEN 61 IF (zh2(i, l2+1) < zhm) extend=.TRUE. 62 END IF 63 IF(extend) THEN 64 l2 = l2+1 65 zsm = zsm + sdsig(l2) 66 zhm = zhm + sdsig(l2) * (zh2(i, l2) - zhm) / zsm 67 END IF 68 END DO 69 70 ! at this point the profile l1-1 (l1 ...l2) l2+1 is stable 71 72 IF(l1<l2) THEN 73 WRITELOG(*,*) 'Mixing between layers ',l1,l2 74 LOG_DBG('convadj') 75 ! perform the mixing of layers l1...l2 : 76 ! * potential temperature is fully mixed, ie replaced by mass-weighted average 77 ! * momentum u,v is mixed with weight zalpha, i.e. u:=zalpha*mean(u)+(1-zalpha)*u 78 ! where zalpha = sum( abs(theta-mean(theta))*dsig) / mean(theta)*sum(dsig) 79 ! for large deviations of theta from its mean, this formula could in principe yield zalpha>1. 80 zalpha=0. 81 zum=0. 82 zvm=0. 83 DO l = l1, l2 84 zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l) 85 zh2(i, l) = zhm 86 ! we must use compute zum, zvm from zu2, zv2 to conserve momentum (see below) 87 zum=zum+dsig(l)*zu2(i,l) 88 zvm=zvm+dsig(l)*zv2(i,l) 89 END DO 90 zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1))) 91 zum=zum/(sig(l1)-sig(l2+1)) 92 zvm=zvm/(sig(l1)-sig(l2+1)) 93 IF(zalpha.GT.1.) THEN 94 WRITELOG(*,*) 'zalpha=',zalpha 95 CALL log_gridpoint(i) 96 LOG_WARN('convadj') 97 STOP 98 END IF 99 IF(zalpha.LT.1.e-5) zalpha=1.e-5 100 ! zum, zvm are mass-weighted averages of zu2, zv2 => mixing preserves total momentum 101 DO l=l1,l2 102 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l)) 103 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) 104 ENDDO 105 END IF 106 107 END DO 108 109 END SUBROUTINE adjust_onecolumn 5 110 6 111 SUBROUTINE convadj(ngrid,nlay,ptimestep, & … … 9 114 & pdufi,pdvfi,pdhfi, & 10 115 & pduadj,pdvadj,pdhadj) 11 USE phys_const12 116 13 117 !======================================================================= 14 118 ! 15 ! ajustement convectif sec 16 ! on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement 119 ! dry convective adjustment 120 ! h = pot. temp. , static instability if ph(above)<ph(below) 121 ! tendencies pdfX are first added to profiles pX, X=u,v,h, yieldig pX2 122 ! adjustment is performed on profiles pX2 123 ! then the difference between adjusted and non-adujsted profiles 124 ! is converted back as tendencies 17 125 ! 18 126 !======================================================================= 19 127 20 !----------------------------------------------------------------------- 21 ! declarations: 22 ! ------------- 128 INTEGER, INTENT(IN) :: ngrid,nlay 129 REAL, INTENT(IN) :: ptimestep 130 REAL, INTENT(IN) :: pplay(ngrid,nlay), pplev(ngrid,nlay+1), ppopsk(ngrid,nlay) 131 REAL, INTENT(IN) :: ph(ngrid,nlay), pu(ngrid,nlay), pv(ngrid,nlay), & 132 & pdhfi(ngrid,nlay), pdufi(ngrid,nlay), pdvfi(ngrid,nlay) 133 REAL, INTENT(OUT) :: pdhadj(ngrid,nlay), pduadj(ngrid,nlay), pdvadj(ngrid,nlay) 134 135 ! local: 136 REAL sig(nlay+1),sdsig(nlay),dsig(nlay) 137 REAL zu(ngrid,nlay), zv(ngrid,nlay), zh(ngrid,nlay) 138 REAL zu2(ngrid,nlay), zv2(ngrid,nlay), zh2(ngrid,nlay) 139 LOGICAL vtest(ngrid) 140 INTEGER jadrs(ngrid) 141 INTEGER jcnt, ig, l, jj 23 142 24 #include "dimensions.h" 143 ! Add physics tendencies to u,v,h 144 DO l=1,nlay 145 DO ig=1,ngrid 146 zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep 147 zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep 148 zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep 149 END DO 150 END DO 151 zu2(:,:)=zu(:,:) 152 zv2(:,:)=zv(:,:) 153 zh2(:,:)=zh(:,:) 25 154 26 ! arguments:27 ! ----------28 29 INTEGER ngrid,nlay30 REAL ptimestep31 REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay)32 REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay)33 REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay)34 REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay)35 36 ! local:37 ! ------38 39 INTEGER ig,i,l,l1,l2,jj40 INTEGER jcnt, jadrs(ngrid)41 42 REAL*8 sig(nlay+1),sdsig(nlay),dsig(nlay)43 REAL*8 zu(ngrid,nlay),zv(ngrid,nlay)44 REAL*8 zh(ngrid,nlay)45 REAL*8 zu2(ngrid,nlay),zv2(ngrid,nlay)46 REAL*8 zh2(ngrid,nlay)47 REAL*8 zhm,zsm,zum,zvm,zalpha48 49 LOGICAL vtest(ngrid),down50 51 !52 !-----------------------------------------------------------------------53 ! initialisation:54 ! ---------------55 !56 !57 155 !----------------------------------------------------------------------- 58 156 ! detection des profils a modifier: … … 63 161 ! sinon, il reste a sa valeur initiale (.FALSE.) 64 162 ! cette operation est vectorisable 65 ! On en profite pour copier la valeur initiale de "ph"66 ! dans le champ de travail "zh"67 163 68 DO l=1,nlay 69 DO ig=1,ngrid 70 zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep 71 zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep 72 zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep 73 END DO 74 END DO 75 76 zu2(:,:)=zu(:,:) 77 zv2(:,:)=zv(:,:) 78 zh2(:,:)=zh(:,:) 79 80 DO ig=1,ngrid 81 vtest(ig)=.FALSE. 82 END DO 83 ! 164 vtest(:)=.FALSE. 84 165 DO l=2,nlay 85 166 DO ig=1,ngrid 86 !CRAY vtest(ig)=CVMGM(.TRUE. , vtest(ig),87 !CRAY . zh2(ig,l)-zh2(ig,l-1))88 167 IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE. 89 168 END DO 90 169 END DO 91 ! 92 !CRAY CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt) 93 jcnt=0 170 171 !---------------------------------- 172 ! Ajustement des profils instables 173 ! -------------------------------- 174 94 175 DO ig=1,ngrid 95 176 IF(vtest(ig)) THEN 96 jcnt=jcnt+197 jadrs(jcnt)=ig177 CALL sigma_levels(ngrid, nlay, ig, pplev, ppopsk, sig, dsig, sdsig) 178 CALL adjust_onecolumn(ngrid, nlay, ig, sig, dsig, sdsig, zu2, zv2, zh2) 98 179 ENDIF 99 180 END DO 100 101 !----------------------------------------------------------------------- 102 ! Ajustement des "jcnt" profils instables indices par "jadrs": 103 ! ------------------------------------------------------------ 104 ! 105 DO jj = 1, jcnt 106 ! 107 i = jadrs(jj) 108 ! 109 ! Calcul des niveaux sigma sur cette colonne 110 DO l=1,nlay+1 111 sig(l)=pplev(i,l)/pplev(i,1) 112 ENDDO 113 DO l=1,nlay 114 dsig(l)=sig(l)-sig(l+1) 115 sdsig(l)=ppopsk(i,l)*dsig(l) 116 ENDDO 117 l2 = 1 118 ! 119 ! -- boucle de sondage vers le haut 120 ! 121 DO WHILE(.TRUE.) 122 ! 8000 CONTINUE 123 ! 124 l2 = l2 + 1 125 ! 126 IF (l2 .GT. nlay) EXIT ! Goto 8001 127 ! 128 IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN 129 ! 130 ! -- l2 est le niveau le plus haut de la colonne instable 131 ! 132 l1 = l2 - 1 133 l = l1 134 zsm = sdsig(l2) 135 zhm = zh2(i, l2) 136 ! 137 ! -- boucle de sondage vers le bas 138 ! 139 DO WHILE(.TRUE.) 140 ! 8020 CONTINUE 141 zsm = zsm + sdsig(l) 142 zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm 143 ! 144 ! -- doit on etendre la colonne vers le bas ? 145 ! 146 !_EC (M1875) 20/6/87 : AND -> AND THEN 147 ! 148 down = .FALSE. 149 IF (l1 .NE. 1) THEN !-- and then 150 IF (zhm .LT. zh2(i, l1-1)) THEN 151 down = .TRUE. 152 END IF 153 END IF 154 ! 155 IF (down) THEN 156 l1 = l1 - 1 157 l = l1 158 ELSE 159 ! -- peut on etendre la colonne vers le haut ? 160 IF (l2 .EQ. nlay) EXIT !Goto 8021 161 IF (zh2(i, l2+1) .GE. zhm) EXIT !Goto 8021 162 l2 = l2 + 1 163 l = l2 164 END IF 165 166 ! GO TO 8020 167 END DO 168 ! 8021 CONTINUE 169 ! 170 ! -- nouveau profil : constant (valeur moyenne) 171 ! 172 zalpha=0. 173 zum=0. 174 zvm=0. 175 DO l = l1, l2 176 zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l) 177 zh2(i, l) = zhm 178 zum=zum+dsig(l)*zu(i,l) 179 zvm=zvm+dsig(l)*zv(i,l) 180 END DO 181 zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1))) 182 zum=zum/(sig(l1)-sig(l2+1)) 183 zvm=zvm/(sig(l1)-sig(l2+1)) 184 IF(zalpha.GT.1.) THEN 185 PRINT*,'WARNING dans convadj zalpha=',zalpha 186 if(ig.eq.1) then 187 print*,'Au pole nord' 188 elseif (ig.eq.ngrid) then 189 print*,'Au pole sud' 190 else 191 print*,'Point i=', & 192 ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1 193 endif 194 STOP 195 zalpha=1. 196 ELSE 197 ! IF(zalpha.LT.0.) STOP'zalpha=0' 198 IF(zalpha.LT.1.e-5) zalpha=1.e-5 199 ENDIF 200 DO l=l1,l2 201 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l)) 202 zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) 203 ENDDO 204 205 l2 = l2 + 1 206 ! 207 END IF 208 ! 209 ! GO TO 8000 181 182 DO l=1,nlay 183 DO ig=1,ngrid 184 pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep 185 pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep 186 pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep 210 187 END DO 211 ! 8001 CONTINUE212 !213 !214 DO l=1,nlay215 DO ig=1,ngrid216 pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep217 pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep218 pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep219 ! pdhadj(ig,l)=0.220 ! pduadj(ig,l)=0.221 ! pdvadj(ig,l)=0.222 END DO223 END DO224 !225 188 END DO 226 189 -
dynamico_lmdz/simple_physics/phyparam/physics/logging.F90
r4199 r4206 1 1 MODULE logging 2 2 3 ! see also logging.h3 ! see also use_logging.h 4 4 ! macro LOGBUF accumulates log output into logging_buffer 5 ! flush_plugin points to a procedure pointer that typically prints all lines in the buffer and prepends tags6 ! This module provides a default implementation of flush_plugin but the top-level driver is welcome to override it.7 5 8 6 IMPLICIT NONE … … 10 8 PRIVATE 11 9 12 INTERFACE 10 INTERFACE ! Explicit interfaces for plugins 11 12 ! Plugin that typically prints all lines in the loggin buffer 'buf' and prepends tags (log level, timestamp, ...) 13 13 SUBROUTINE plugin(lev, tag, buf) 14 14 INTEGER, INTENT(IN) :: lev 15 15 CHARACTER(*), INTENT(IN) :: tag, buf(:) 16 16 END SUBROUTINE plugin 17 18 ! Plugin that writes into string 'line' information about the gridpoint of index 'index' 19 SUBROUTINE plugin_log_gridpoint(index, line) 20 INTEGER, INTENT(IN) :: index ! index of gridpoint 21 CHARACTER(*), INTENT(OUT) :: line 22 END SUBROUTINE plugin_log_gridpoint 23 17 24 END INTERFACE 18 25 19 PROCEDURE(plugin), POINTER :: flush_plugin => default_flush_plugin 26 ! This module provides a default implementation of flush_plugin but the top-level driver is welcome to override it. 27 PROCEDURE(plugin), POINTER :: flush_plugin => default_flush 28 29 ! The top-level driver MUST provide an implementation for log_gridpoint_plugin 30 PROCEDURE(plugin_log_gridpoint), POINTER :: log_gridpoint_plugin => NULL() 20 31 21 32 INTEGER, PARAMETER :: bufsize=10000 … … 24 35 INTEGER :: logging_lineno=0 25 36 26 INTEGER, PARAMETER, PUBLIC :: log_level_ dbg=0, log_level_info=137 INTEGER, PARAMETER, PUBLIC :: log_level_fatal=0, log_level_error=1, log_level_warn=2, log_level_info=3, log_level_dbg=4 27 38 28 PUBLIC :: logging_buf, logging_lineno, flush_log, flush_plugin 39 PUBLIC :: logging_buf, logging_lineno, flush_log, log_gridpoint, & 40 flush_plugin, log_gridpoint_plugin 29 41 30 42 CONTAINS … … 37 49 END SUBROUTINE flush_log 38 50 39 SUBROUTINE default_flush _plugin(lev, tag, buf)51 SUBROUTINE default_flush(lev, tag, buf) 40 52 INTEGER, INTENT(IN) :: lev 41 53 CHARACTER(*), INTENT(IN) :: tag, buf(:) … … 44 56 PRINT *, '[INFO ',tag,']', TRIM(buf(i)) 45 57 END DO 46 END SUBROUTINE default_flush_plugin 47 58 END SUBROUTINE default_flush 59 60 SUBROUTINE log_gridpoint(index) 61 INTEGER, INTENT(IN) :: index 62 logging_lineno = logging_lineno+1 63 CALL log_gridpoint_plugin(index, logging_buf(logging_lineno)) 64 END SUBROUTINE log_gridpoint 65 48 66 END MODULE logging -
dynamico_lmdz/simple_physics/phyparam/physics/use_logging.h
r4199 r4206 1 1 USE logging 2 2 #define WRITELOG(junk, fmt) logging_lineno = logging_lineno+1 ; WRITE(logging_buf(logging_lineno), fmt) 3 #define LOG_WARN(tag) CALL flush_log(log_level_warn, tag) 3 4 #define LOG_INFO(tag) CALL flush_log(log_level_info, tag) 4 5 #define LOG_DBG(tag) CALL flush_log(log_level_dbg, tag)
Note: See TracChangeset
for help on using the changeset viewer.
