Changeset 1260 for trunk/LMDZ.MARS/libf/aeronomars

Timestamp:

May 14, 2014, 3:12:01 PM (11 years ago)

Author:

emillour

Message:

Mars GCM:

• Bug fix in jthermcalc.F: rounding issues in computation of the solar zenith angle when SZA slightly greater than 90
• Corrected misleading warning in nlte_tcool.F

FGG

File:

• trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F (modified) (8 diffs)

trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F

@@ -154 +154 @@
      $        h2colx(i)*crscabsi2(5,indexint) + 
      $        ncolx(i)*crscabsi2(9,indexint)
+         
+         
       end do
       limdown=1.e-20
       limup=1.e26
 
-      ! Ehouarn: sanity check
-      ! test that auxcolinp in monotonously increasing from 1 to nlayermx
-      do j=1,nlayermx-1
-        if (auxcolinp(j).gt.auxcolinp(j+1)) then
-          !there is a problem
-          write(*,*) "jthermcalc error: "
-          write(*,*) "auxcolinp() not increasing with altitude index!"
-          write(*,*) "j=",j," auxcolinp(j)=",auxcolinp(j)
-          write(*,*) "                auxcolinp(j+1)=",auxcolinp(j+1)
-          ! Quick fix:
-          if (j==1) then
-            auxcolinp(j)=auxcolinp(j+1)/2.
-          else
-            ! compute it as a geometric mean from encompassing values
-            auxcolinp(j)=sqrt(auxcolinp(j-1)*auxcolinp(j+1))
-          endif
-          write(*,*) " Quick fixed to auxcolinp(j)=",auxcolinp(j)
-        endif
-      enddo
 
 c     Interpolations
@@ -206 +189 @@
          ind=auxind(i)
          auxi=nlayermx-i+1
+         ! Ehouarn: test
+          if ((ind+1.gt.nz2).or.(ind.le.0)) then
+            write(*,*) "jthercalc error: ind=",ind,ig,zenit
+            write(*,*) " auxind(1:nlayermx)=",auxind
+            write(*,*) " auxcolinp(:nlayermx)=",auxcolinp
+            write(*,*) " co2colx(:nlayermx)=",co2colx
+            write(*,*) " o2colx(:nlayermx)=",o2colx
+            write(*,*) " o3pcolx(:nlayermx)=",o3pcolx
+            write(*,*) " h2colx(:nlayermx)=",h2colx
+            write(*,*) " ncolx(:nlayermx)=",ncolx
+            write(*,*) " auxcoltab(1:nz2)=",auxcoltab
+            write(*,*) " limdown=",limdown
+            write(*,*) " limup=",limup
+            call abort_gcm('jthermcalc','error',1)
+          endif
          !CO2 interpolated coefficient
          jfotsout(indexint,1,auxi) = wm(i)*auxjco2(ind+1) +
@@ -1347 +1345 @@
             end do    !Of do j=1,nlayesp
          endif        !Of ilayesp(nlayesp).eq.-1
+
       enddo           !Of do i=nlayermx,1,-1
+
+
       return
 
@@ -1438 +1439 @@
         real*8      rkmmini               ! distance zmini to center of P [km] 
         real*8      rkmj                  ! intermediate distance to C of P [km]
-
 c external function
         external  grid_R8          ! Returns index of layer containing the altitude
@@ -1536 +1536 @@
 
            zmini=(radio+zz)*sin(szarad)-radio
+           !zmini should be lower than zz, as SZA<90. However, in situations
+           !where SZA is very close to 90, rounding errors can make zmini
+           !slightly higher than zz, causing problems in the determination
+           !of the jmin index. A correction is implemented in the determination
+           !of jmin, some lines below
            rkmmini = radio + zmini
 
@@ -1545 +1550 @@
            else
               jmin=grid_R8(zmini,diz,nlayermx)+1
-              
+              !Correction for possible rounding errors when SZA very close 
+              !to 90 degrees
+              if(jmin.gt.grid_R8(zz,diz,nlayermx)) then
+                 write(*,*)'jthermcalc warning: possible rounding error'
+                 write(*,*)'point,sza,layer:',ig,szadeg,capa
+                 jmin=grid_R8(zz,diz,nlayermx)
+              endif
 
               if(abs(zz-diz(nlayermx)).lt.1.d-4) goto 9876
@@ -1620 +1631 @@
         end if
 
+
         return
 
@@ -1652 +1664 @@
 c*** CODE START 
 
-        if ( z .lt. zgrid(1) .or. z.gt.zgrid(nz) ) then 
+        if ( z .lt. zgrid(1)  ) then 
            write (*,*) ' GRID/ z outside bounds of zgrid '
            write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
-           stop ' Serious error in GRID.F '
+           z = zgrid(1) 
+           write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)'
+           write(*,*) 'Please check values of z and zgrid above'
+        endif
+        if (z .gt. zgrid(nz) ) then
+           write (*,*) ' GRID/ z outside bounds of zgrid '
+           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
+           z = zgrid(nz) 
+           write(*,*) 'WARNING: error in grid_r8 (jthermcalc.F)'
+           write(*,*) 'Please check values of z and zgrid above'
         endif
         if ( nz .lt. 2 ) then