Ignore:
Timestamp:
Dec 20, 2013, 4:04:56 PM (11 years ago)
Author:
emillour
Message:

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

  • Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
  • dyn3d/control.h is now module control_mod.F90
  • rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
  • added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
  • created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

EM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/surfini.F

    r1087 r1130  
    88     &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
    99     &                     albedice
     10      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
     11      use mod_phys_lmdz_para, only : is_master, gather, scatter
    1012      IMPLICIT NONE
    1113c=======================================================================
    … …  
    2729#include "datafile.h"
    2830
    29       INTEGER ngrid,ig,icap,iq,alternate
    30       REAL  piceco2(ngrid),psolaralb(ngrid,2)
    31       REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
     31      integer,intent(in) :: ngrid ! number of atmospheric columns
     32      real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
     33      real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
     34      real,intent(out) :: psolaralb(ngrid,2) ! albedo
     35
     36      INTEGER ig,icap,iq,alternate
    3237      REAL icedryness ! ice dryness
    3338     
    3439      ! longwatercaptag is watercaptag. Trick for some compilers
    3540      LOGICAL, DIMENSION(100000) :: longwatercaptag
    36 
    37       EXTERNAL ISMIN,ISMAX
    38       INTEGER ISMIN,ISMAX
    3941     
    4042! There are 3 different modes for ice distribution:
    … …  
    5355      REAL        zelat,zelon
    5456
     57#ifdef CPP_PARA
    5558      INTEGER nb_ice(ngrid,2)              ! number of counts | detected ice for GCM grid
     59#else
     60      INTEGER nb_ice(klon_glo,2)
     61#endif
    5662      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
    5763
    … …  
    6470     
    6571      character (len=100) :: zedatafile
     72
     73      ! to handle parallel cases
     74#if CPP_PARA
     75      logical watercaptag_glo(klon_glo)
     76      real dryness_glo(klon_glo)
     77      real lati_glo(klon_glo)
     78      real long_glo(klon_glo)
     79#else
     80      logical watercaptag_glo(ngrid)
     81      real dryness_glo(ngrid)
     82      real lati_glo(ngrid)
     83      real long_glo(ngrid)
     84#endif
     85
    6686c
    6787c=======================================================================
     88! Initialize watercaptag (default is false)
     89      watercaptag_glo(:)=.false.
    6890
    6991c     water ice outliers
    … …  
    93115         endif
    94116         
    95          write(*,*) "Ice dryness ?"
     117         write(*,*) "surfini: Ice dryness ?"
    96118         icedryness=1. ! default value
    97119         call getin("icedryness",icedryness)
    98          write(*,*) " icedryness = ",icedryness
     120         write(*,*) "surfini: icedryness = ",icedryness
    99121         dryness (:) = icedryness
    100122         
    … …  
    125147#else
    126148
    127 
    128 
    129       IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
     149      ! To be able to run in parallel, we work on the full grid
     150      ! and dispatch results afterwards
     151
     152      ! start by geting latitudes and logitudes on full grid
     153      ! (in serial mode, this is just a copy)
     154      call gather(lati,lati_glo)
     155      call gather(long,long_glo)
     156
     157      if (is_master) then
     158
     159        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
    130160     
    131161         print*, 'ngrid = 1, do no put ice caps in surfini.F'
    132162
    133       ELSE IF (icelocationmode .eq. 1) THEN
     163        ELSE IF (icelocationmode .eq. 1) THEN
    134164     
    135165         print*,'Surfini: ice caps defined from surface.nc'
    … …  
    144174
    145175           
    146         zedatafile = trim(datafile)
     176         zedatafile = trim(datafile)
    147177 
    148178       
    149         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
    150      &  NF90_NOWRITE,nid)
     179         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
     180     &   NF90_NOWRITE,nid)
    151181     
    152       IF (ierr.NE.nf90_noerr) THEN
     182         IF (ierr.NE.nf90_noerr) THEN
    153183       write(*,*)'Error : cannot open file surface.nc '
    154184       write(*,*)'(in phymars/surfini.F)'
    … …  
    160190       write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
    161191       CALL ABORT
    162       ENDIF
    163      
    164      
    165       ierr=nf90_inq_varid(nid, string, nvarid)
    166       if (ierr.ne.nf90_noerr) then
    167         write(*,*) 'surfini error, cannot find ',trim(string)
    168         write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
    169         write(*,*)trim(nf90_strerror(ierr))
    170         stop
    171       endif
    172 
    173       ierr=nf90_get_var(nid, nvarid, zdata)
    174 
    175       if (ierr.ne.nf90_noerr) then
    176         write(*,*) 'surfini: error failed loading ',trim(string)
    177         write(*,*)trim(nf90_strerror(ierr))
    178         stop
    179       endif
     192         ENDIF
     193     
     194     
     195         ierr=nf90_inq_varid(nid, string, nvarid)
     196         if (ierr.ne.nf90_noerr) then
     197          write(*,*) 'surfini error, cannot find ',trim(string)
     198          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
     199          write(*,*)trim(nf90_strerror(ierr))
     200          stop
     201         endif
     202
     203         ierr=nf90_get_var(nid, nvarid, zdata)
     204
     205         if (ierr.ne.nf90_noerr) then
     206          write(*,*) 'surfini: error failed loading ',trim(string)
     207          write(*,*)trim(nf90_strerror(ierr))
     208          stop
     209         endif
    180210 
    181211                     
    182       ierr=nf90_close(nid)
     212         ierr=nf90_close(nid)
    183213 
    184214
    185       nb_ice(:,1) = 1 ! default: there is no ice
    186       latice(:,1) = 1
    187       lonice(:,1) = 1
    188       nb_ice(:,2) = 0
    189       latice(:,2) = 0
    190       lonice(:,2) = 0
    191       !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
    192 
    193       ! loop over the GCM grid - except for poles (ig=1 and ngrid)
    194       do ig=2,ngrid-1
    195      
    196         ! loop over the surface file grid     
    197         do i=1,imd
    198           do j=1,jmd
    199             zelon = i - 180.
    200             zelat = 90. - j
    201             if ((abs(lati(ig)*180./pi-zelat) .le. 90./real(jjm)) .and.
    202      &        (abs(long(ig)*180./pi-zelon) .le. 180./real(iim))) then
     215         nb_ice(:,1) = 1 ! default: there is no ice
     216         latice(:,1) = 1
     217         lonice(:,1) = 1
     218         nb_ice(:,2) = 0
     219         latice(:,2) = 0
     220         lonice(:,2) = 0
     221         !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
     222
     223         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
     224         do ig=2,klon_glo-1
     225     
     226          ! loop over the surface file grid     
     227          do i=1,imd
     228           do j=1,jmd
     229             zelon = i - 180.
     230             zelat = 90. - j
     231            if ((abs(lati_glo(ig)*180./pi-zelat).le.90./real(jjm)) .and.
     232     &        (abs(long_glo(ig)*180./pi-zelon).le.180./real(iim))) then
    203233              ! count all points in that GCM grid point
    204234              nb_ice(ig,1) = nb_ice(ig,1) + 1
    … …  
    207237     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
    208238             endif
    209           enddo
    210         enddo 
     239           enddo
     240          enddo 
    211241
    212242        ! projection of nb_ice on GCM lat and lon axes
    213         latice(1+(ig-2)/iim,:) =
     243          latice(1+(ig-2)/iim,:) =
    214244     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
    215         lonice(1+mod(ig-2,iim),:) =
     245          lonice(1+mod(ig-2,iim),:) =
    216246     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
    217247
    218       enddo ! of do ig=2,ngrid-1
     248         enddo ! of do ig=2,klon_glo-1
    219249     
    220250
    221251     
    222       ! special case for poles
    223       nb_ice(1,2)   = 1  ! ice prescribed on north pole
    224       latice(1,:)   = nb_ice(1,:)
    225       lonice(1,:)   = nb_ice(1,:)
    226       latice(jjm,:) = nb_ice(ngrid,:)
    227       lonice(iim,:) = nb_ice(ngrid,:)
     252         ! special case for poles
     253         nb_ice(1,2)   = 1  ! ice prescribed on north pole
     254         latice(1,:)   = nb_ice(1,:)
     255         lonice(1,:)   = nb_ice(1,:)
     256         latice(jjm,:) = nb_ice(ngrid,:)
     257         lonice(iim,:) = nb_ice(ngrid,:)
    228258     
    229259     
    … …  
    241271     
    242272   
    243       ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
    244       do i=1,jjm/2
    245       step  = 1. ! threshold to add ice cap
    246       count = 0. ! number of ice GCM caps at this latitude
    247       ! ratiolat is the ratio of area covered by ice within this GCM latitude range
    248       ratiolat  = real(latice(i,2))/real(latice(i,1))
    249       !print*,'i',i,(i-1)*iim+2,i*iim+1
     273         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
     274         do i=1,jjm/2
     275          step  = 1. ! threshold to add ice cap
     276          count = 0. ! number of ice GCM caps at this latitude
     277          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
     278          ratiolat  = real(latice(i,2))/real(latice(i,1))
     279          !print*,'i',i,(i-1)*iim+2,i*iim+1
    250280     
    251         ! put ice caps while there is not enough ice,
    252         ! as long as the threshold is above 20%
    253         do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
    254           count = 0.
    255           ! loop over GCM longitudes
    256           do j=1,iim
     281          ! put ice caps while there is not enough ice,
     282          ! as long as the threshold is above 20%
     283          do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
     284           count = 0.
     285           ! loop over GCM longitudes
     286           do j=1,iim
    257287            ! if the detected ice ratio in the GCM grid point
    258288            ! is more than 'step', then add ice
    259289            if (real(nb_ice((i-1)*iim+1+j,2))
    260290     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
    261                   watercaptag((i-1)*iim+1+j) = .true.
     291                  watercaptag_glo((i-1)*iim+1+j) = .true.
    262292                  count = count + 1
    263293            endif
    264           enddo ! of do j=1,iim
     294           enddo ! of do j=1,iim
     295           !print*, 'step',step,count,ratiolat*iim
     296           step = step - 0.01
     297          enddo ! of do while
    265298          !print*, 'step',step,count,ratiolat*iim
    266           step = step - 0.01
    267         enddo ! of do while
    268       !print*, 'step',step,count,ratiolat*iim
    269 
    270       enddo ! of do i=1,jjm/2
     299
     300         enddo ! of do i=1,jjm/2
    271301           
    272302
    273       ELSE IF (icelocationmode .eq. 2) THEN
    274      
    275         print*,'Surfini: predefined ice caps'
    276      
    277         if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
     303        ELSE IF (icelocationmode .eq. 2) THEN
     304     
     305         print*,'Surfini: predefined ice caps'
     306     
     307         if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
    278308           
    279309          print*,'water ice caps distribution for 32x24 resolution'
    … …  
    288318!---------------------   OUTLIERS  ----------------------------
    289319
    290         else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
     320         else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
    291321
    292322          print*,'water ice caps distribution for 64x48 resolution'
    … …  
    323353
    324354           
    325         else if (ngrid .ne. 1) then
     355         else if (klon_glo .ne. 1) then
    326356       
    327       print*,'No predefined ice location for this resolution :',iim,jjm
    328       print*,'Please change icelocationmode in surfini.F'
    329       print*,'Or add some new definitions ...'
    330       call abort
     357          print*,'No predefined ice location for this resolution :',
     358     &           iim,jjm
     359          print*,'Please change icelocationmode in surfini.F'
     360          print*,'Or add some new definitions ...'
     361          call abort
    331362         
    332         endif
    333 
    334         do ig=1,ngrid
    335           if (longwatercaptag(ig)) watercaptag(ig) = .true.
    336         enddo
    337 
    338 
    339       ELSE IF (icelocationmode .eq. 3) THEN
    340      
    341         print*,'Surfini: ice caps defined by lat and lon values'
    342 
    343          do ig=1,ngrid
     363         endif
     364
     365         do ig=1,klon_glo
     366          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
     367         enddo
     368
     369
     370        ELSE IF (icelocationmode .eq. 3) THEN
     371     
     372         print*,'Surfini: ice caps defined by lat and lon values'
     373
     374         do ig=1,klon_glo
    344375         
    345376c-------- Towards olympia planitia water caps -----------
    346377c--------------------------------------------------------
    347378
    348        if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
    349      .           ( lati(ig)*180./pi .le. 80.  ) .and.
    350      .           ( long(ig)*180./pi .ge. 110. ) .and.
    351      .           ( long(ig)*180./pi .le. 181. ) )
     379          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
     380     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
     381     .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
     382     .           ( long_glo(ig)*180./pi .le. 181. ) )
    352383     .         .or.
    353384
    354      .         ( ( lati(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
    355      .           ( lati(ig)*180./pi .le. 76.  ) .and.
    356      .           ( long(ig)*180./pi .ge. 150. ) .and.
    357      .           ( long(ig)*180./pi .le. 168. ) )
     385     .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
     386     .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
     387     .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
     388     .           ( long_glo(ig)*180./pi .le. 168. ) )
    358389     .         .or.
    359      .         ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
    360      .           ( lati(ig)*180./pi .le. 80.  ) .and.
    361      .           ( long(ig)*180./pi .ge. -150.) .and.
    362      .           ( long(ig)*180./pi .le. -110.) ) )
     390     .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
     391     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
     392     .           ( long_glo(ig)*180./pi .ge. -150.) .and.
     393     .           ( long_glo(ig)*180./pi .le. -110.) ) )
    363394     .         then
    364395             
    … …  
    370401               endif !end if alternate = 0
    371402               
    372        endif
     403          endif
    373404
    374405c----------- Opposite olympia planitia water cap --------
    375406c--------------------------------------------------------
    376407
    377         if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
    378      .         ( lati(ig)*180./pi     .le.  84 ) )
     408          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
     409     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
    379410     .         .and.
    380      .       ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
    381      .         ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
    382 !     .     ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
    383 !     .         ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
    384 !     .       ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
    385 !     .         ( long(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
    386         !   watercaptag(ig)=.true.
    387         endif
     411     .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
     412     .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
     413!     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
     414!     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
     415!     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
     416!     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
     417        !   watercaptag_glo(ig)=.true.
     418          endif
    388419
    389420
    … …  
    391422c--------------------------------------------------------
    392423
    393       if (abs(lati(ig)*180./pi).gt.80)
    394      .       watercaptag(ig)=.true.
     424          if (abs(lati_glo(ig)*180./pi).gt.80)
     425     .          watercaptag_glo(ig)=.true.
    395426           
    396427c--------------------------------------------------------
    397428c--------------------------------------------------------
    398       end do ! of (ngrid)
    399 
    400 
    401        ELSE
     429         end do ! of (klon_glo)
     430
     431
     432        ELSE
    402433     
    403434         print*, 'In surfini.F, icelocationmode is ', icelocationmode
    … …  
    405436         call abort
    406437
    407        ENDIF ! of if (icelocation)
     438        ENDIF ! of if (icelocation)
    408439       
    409440       
    410441        ! print caps locations - useful for plots too
    411         print*,'latitude | longitude | ig'
    412         do ig=1,ngrid
    413           dryness (ig) = icedryness
    414 
    415           if (watercaptag(ig)) then
    416              print*,'ice water cap', lati(ig)*180./pi,
    417      .              long(ig)*180./pi, ig
     442        print*,'surfini: latitude | longitude | ig'
     443        do ig=1,klon_glo
     444          dryness_glo(ig) = icedryness
     445
     446          if (watercaptag_glo(ig)) then
     447             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
     448     &              long_glo(ig)*180./pi, ig
    418449          endif
    419450        enddo
    420451       
     452       endif !of if (is_master)
     453       
     454       ! Now scatter fields watercaptag and dryness from master to all
     455       ! (is just a plain copy in serial mode)
     456       call scatter(dryness_glo,dryness)
     457       call scatter(watercaptag_glo,watercaptag)
     458       
    421459#endif     
    422 
     460! end of #else of #ifdef MESOSCALE
    423461       ENDIF ! (caps & water)
    424462       
    … …  
    439477              psolaralb(ig,2)=albedodat(ig)
    440478         END DO
    441          PRINT*,'minimum albedo sans water caps',
    442      s     albedodat(ISMIN(ngrid,albedodat,1))
    443          PRINT*,'maximum albedo sans water caps',
    444      s     albedodat(ISMAX(ngrid,albedodat,1))
     479         PRINT*,'surfini: minimum albedo without water caps',
     480     &          minval(albedodat)
     481         PRINT*,'surfini: maximum albedo without water caps',
     482     &          maxval(albedodat)
    445483
    446484c        initial albedo if permanent H2O ice is present
    … …  
    453491            ENDIF
    454492           END DO
    455            PRINT*,'minimum albedo avec water caps',
    456      s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
    457            PRINT*,'maximum albedo avec water caps',
    458      s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
     493           PRINT*,'surfini: minimum albedo with water caps',
     494     &            minval(albedodat)
     495           PRINT*,'surfini: maximum albedo with water caps',
     496     &            maxval(albedodat)
    459497
    460498         ENDIF
    … …  
    465503       DO ig=1,ngrid
    466504         IF (piceco2(ig) .GT. 0.) THEN
    467              IF(ig.GT.ngrid/2+1) THEN
    468                 icap=2
     505             IF(lati(ig).LT. 0.) THEN
     506                icap=2 ! Southern hemisphere
    469507             ELSE
    470                 icap=1
     508                icap=1 ! Northern hemisphere
    471509             ENDIF
    472510             psolaralb(ig,1) = albedice(icap)
    … …  
    514552           endif
    515553          end do
    516           PRINT*,'minimum albedo avec givre et co2',
    517      s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
    518           PRINT*,'maximum albedo avec givre et co2',
    519      s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
     554          PRINT*,'surfini: minimum albedo with frost and co2',
     555     &            minval(albedodat)
     556          PRINT*,'surfini: maximum albedo with frost and co2',
     557     &            maxval(albedodat)
    520558       ELSE
    521559         watercaptag(:) = .false.
    522        END IF
     560       END IF ! OF IF(water)
    523561         
    524562      RETURN
Note: See TracChangeset for help on using the changeset viewer.