source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 11 years ago

Mars GCM: (a first step towards using parallel dynamics)

  • IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
    • advtrac.h (in dyn3d) removed and replaced by module infotrac.F
    • tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
  • Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
  • Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

EM
File size: 18.1 KB
Line 
1      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
2   ! to use  'getin'
3      USE ioipsl_getincom
4      use netcdf
5      use tracer_mod, only: nqmx, noms, dryness
6      IMPLICIT NONE
7c=======================================================================
8c
9c   creation des calottes pour l'etat initial
10c
11c=======================================================================
12c-----------------------------------------------------------------------
13c   Declarations:
14c   -------------
15#include "dimensions.h"
16#include "dimphys.h"
17#include "surfdat.h"
18#include "callkeys.h"
19!#include "tracer.h"
20#include "comgeomfi.h"
21#include "comcstfi.h"
22
23#include "datafile.h"
24
25      INTEGER ngrid,ig,icap,iq,alternate
26      REAL  piceco2(ngrid),psolaralb(ngrid,2)
27      REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
28      REAL icedryness ! ice dryness
29     
30      ! longwatercaptag is watercaptag. Trick for some compilers
31      LOGICAL, DIMENSION(100000) :: longwatercaptag
32
33      EXTERNAL ISMIN,ISMAX
34      INTEGER ISMIN,ISMAX
35     
36! There are 3 different modes for ice distribution:
37! icelocationmode = 1 ---> based on data from surface.nc
38! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
39! icelocationmode = 3 ---> based on logical relations for latitude and longitude
40! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
41      INTEGER,SAVE :: icelocationmode = 2
42       
43       
44      !in case icelocationmode == 1
45      INTEGER i,j
46      INTEGER     imd,jmd
47      PARAMETER   (imd=360,jmd=180)
48      REAL        zdata(imd,jmd)
49      REAL        zelat,zelon
50
51      INTEGER nb_ice(ngrid,2)              ! number of counts | detected ice for GCM grid
52      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
53
54      REAL step,count,ratiolat
55
56      INTEGER   ierr,nid,nvarid
57     
58      REAL,SAVE :: min_icevalue = 500.
59      character(len=50) :: string = 'thermal'
60     
61      character (len=100) :: zedatafile
62c
63c=======================================================================
64
65c     water ice outliers
66c     ------------------------------------------
67
68      IF ((water) .and. (caps)) THEN
69     
70c Perennial H20 north cap defined by watercaptag=true (allows surface to be
71c hollowed by sublimation in vdifc).
72
73c We might not want albedodat to be modified because it is used to write
74c restart files. Instead, albedo is directly modified when needed (i.e.
75c if we have watercaptag and no co2 ice), below and in albedocaps.F90
76
77c       "Dryness coefficient" controlling the evaporation and
78c        sublimation from the ground water ice (close to 1)
79c        HERE, the goal is to correct for the fact
80c        that the simulated permanent water ice polar caps
81c        is larger than the actual cap and the atmospheric
82c        opacity not always realistic.
83
84         alternate = 0
85         
86         if (ngridmx .ne. 1) then
87           watercaptag(:) = .false.
88           longwatercaptag(:) = .false.
89         endif
90         
91         write(*,*) "Ice dryness ?"
92         icedryness=1. ! default value
93         call getin("icedryness",icedryness)
94         write(*,*) " icedryness = ",icedryness
95         dryness (:) = icedryness
96         
97       
98#ifdef MESOSCALE
99
100      do ig=1,ngridmx
101
102         !write(*,*) "all qsurf to zero. dirty."
103         do iq=1,nqmx
104         qsurf(ig,iq)=0.  !! on jette les inputs GCM
105                          !! on regle juste watercaptag
106                          !! il faudrait garder les inputs GCM
107                          !! si elles sont consequentes
108         enddo
109         if ( ( lati(ig)*180./pi .gt. 70. ) .and.
110     .        ( albedodat(ig) .ge. 0.26   ) )  then
111                 write(*,*)"outlier ",ig
112                 watercaptag(ig)  = .true.
113                 dryness(ig)      = 1.
114                 albedodat(ig)    = albedo_h2o_ice  !! pour output
115         else
116                 watercaptag(ig)  = .false.
117                 dryness(ig)      = 1.
118         endif
119         
120      enddo
121#else
122
123
124
125      IF (ngridmx .eq. 1) THEN ! special case for 1d --> do nothing
126     
127         print*, 'ngridmx = 1, do no put ice caps in surfini.F'
128
129      ELSE IF (icelocationmode .eq. 1) THEN
130     
131         print*,'Surfini: ice caps defined from surface.nc'
132           
133! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
134! Typically, it is for thermal inertia above 500 tiu.
135! Two conditions are verified:
136! 1. GCM ice caps are defined such as area is conserved for a given latitude
137! (the approximation is that all points within the GCM latitude resolution have the same area).
138! 2. caps are placed to fill the GCM points with the most detected ice first.
139     
140
141           
142        zedatafile = trim(datafile)
143 
144       
145        ierr=nf90_open(trim(zedatafile)//'/surface.nc',
146     &  NF90_NOWRITE,nid)
147     
148      IF (ierr.NE.nf90_noerr) THEN
149       write(*,*)'Error : cannot open file surface.nc '
150       write(*,*)'(in phymars/surfini.F)'
151       write(*,*)'It should be in :',trim(zedatafile),'/'
152       write(*,*)'1) You can set this path in the callphys.def file:'
153       write(*,*)'   datadir=/path/to/the/datafiles'
154       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
155       write(*,*)'   can be obtained online on:'
156       write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
157       CALL ABORT
158      ENDIF
159     
160     
161      ierr=nf90_inq_varid(nid, string, nvarid)
162      if (ierr.ne.nf90_noerr) then
163        write(*,*) 'surfini error, cannot find ',trim(string)
164        write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
165        write(*,*)trim(nf90_strerror(ierr))
166        stop
167      endif
168
169      ierr=nf90_get_var(nid, nvarid, zdata)
170
171      if (ierr.ne.nf90_noerr) then
172        write(*,*) 'surfini: error failed loading ',trim(string)
173        write(*,*)trim(nf90_strerror(ierr))
174        stop
175      endif
176 
177                     
178      ierr=nf90_close(nid)
179 
180
181      nb_ice(:,1) = 1 ! default: there is no ice
182      latice(:,1) = 1
183      lonice(:,1) = 1
184      nb_ice(:,2) = 0
185      latice(:,2) = 0
186      lonice(:,2) = 0
187      !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
188
189      ! loop over the GCM grid - except for poles (ig=1 and ngridmx)
190      do ig=2,ngridmx-1
191     
192        ! loop over the surface file grid     
193        do i=1,imd
194          do j=1,jmd
195            zelon = i - 180.
196            zelat = 90. - j
197            if ((abs(lati(ig)*180./pi-zelat) .le. 90./real(jjm)) .and.
198     &        (abs(long(ig)*180./pi-zelon) .le. 180./real(iim))) then
199              ! count all points in that GCM grid point
200              nb_ice(ig,1) = nb_ice(ig,1) + 1
201              if (zdata(i,j) > min_icevalue)
202                 ! count all detected points in that GCM grid point
203     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
204             endif
205          enddo
206        enddo 
207
208        ! projection of nb_ice on GCM lat and lon axes
209        latice(1+(ig-2)/iim,:) =
210     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
211        lonice(1+mod(ig-2,iim),:) =
212     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
213
214      enddo ! of do ig=2,ngridmx-1
215     
216
217     
218      ! special case for poles
219      nb_ice(1,2)   = 1  ! ice prescribed on north pole
220      latice(1,:)   = nb_ice(1,:)
221      lonice(1,:)   = nb_ice(1,:)
222      latice(jjm,:) = nb_ice(ngridmx,:)
223      lonice(iim,:) = nb_ice(ngridmx,:)
224     
225     
226!      print*, 'latice TOT', latice(:,1)
227!      print*, 'latice FOUND', latice(:,2)
228!      print*, 'lonice TOT', lonice(:,1)
229!      print*, 'lonice FOUND', lonice(:,2)
230     
231!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
232!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
233     
234!      print*,''
235!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
236!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
237     
238   
239      ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
240      do i=1,jjm/2
241      step  = 1. ! threshold to add ice cap
242      count = 0. ! number of ice GCM caps at this latitude
243      ! ratiolat is the ratio of area covered by ice within this GCM latitude range
244      ratiolat  = real(latice(i,2))/real(latice(i,1))
245      !print*,'i',i,(i-1)*iim+2,i*iim+1
246     
247        ! put ice caps while there is not enough ice,
248        ! as long as the threshold is above 20%
249        do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
250          count = 0.
251          ! loop over GCM longitudes
252          do j=1,iim
253            ! if the detected ice ratio in the GCM grid point
254            ! is more than 'step', then add ice
255            if (real(nb_ice((i-1)*iim+1+j,2))
256     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
257                  watercaptag((i-1)*iim+1+j) = .true.
258                  count = count + 1
259            endif
260          enddo ! of do j=1,iim
261          !print*, 'step',step,count,ratiolat*iim
262          step = step - 0.01
263        enddo ! of do while
264      !print*, 'step',step,count,ratiolat*iim
265
266      enddo ! of do i=1,jjm/2
267           
268
269      ELSE IF (icelocationmode .eq. 2) THEN
270     
271        print*,'Surfini: predefined ice caps'
272     
273        if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
274           
275          print*,'water ice caps distribution for 32x24 resolution'
276          longwatercaptag(1:9)    = .true. ! central cap - core
277          longwatercaptag(26:33)  = .true. ! central cap
278          longwatercaptag(1:33)  = .true. ! central cap
279          longwatercaptag(56)  = .true. ! central cap
280          longwatercaptag(58)  = .true. ! central cap
281          longwatercaptag(60)  = .true. ! central cap
282          longwatercaptag(62)  = .true. ! central cap
283          longwatercaptag(64)  = .true. ! central cap
284!---------------------   OUTLIERS  ----------------------------
285
286        else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
287
288          print*,'water ice caps distribution for 64x48 resolution'
289          longwatercaptag(1:65)   = .true. ! central cap - core
290          longwatercaptag(75:85)  = .true. ! central cap
291          longwatercaptag(93:114) = .true. ! central cap
292!---------------------   OUTLIERS  ----------------------------
293          if (.true.) then
294          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
295          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
296          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
297          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
298          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
299          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
300          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
301          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
302          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
303          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
304          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
305          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
306          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
307          longwatercaptag(194)    = .true. ! outlier, lat = 75
308          longwatercaptag(203)    = .true. ! outlier, lat = 75
309          longwatercaptag(207)    = .true. ! outlier, lat = 75
310          longwatercaptag(244)    = .true. ! outlier, lat = 75
311          longwatercaptag(246)    = .true. ! outlier, lat = 75
312          longwatercaptag(250)    = .true. ! outlier, lat = 75
313          longwatercaptag(252)    = .true. ! outlier, lat = 75
314          longwatercaptag(254)    = .true. ! outlier, lat = 75
315          longwatercaptag(256)    = .true. ! outlier, lat = 75
316          endif
317!--------------------------------------------------------------       
318
319
320           
321        else if (ngridmx .ne. 1) then
322       
323      print*,'No predefined ice location for this resolution :',iim,jjm
324      print*,'Please change icelocationmode in surfini.F'
325      print*,'Or add some new definitions ...'
326      call abort
327         
328        endif
329
330        do ig=1,ngridmx
331          if (longwatercaptag(ig)) watercaptag(ig) = .true.
332        enddo
333
334
335      ELSE IF (icelocationmode .eq. 3) THEN
336     
337        print*,'Surfini: ice caps defined by lat and lon values'
338
339         do ig=1,ngridmx
340         
341c-------- Towards olympia planitia water caps -----------
342c--------------------------------------------------------
343
344       if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
345     .           ( lati(ig)*180./pi .le. 80.  ) .and.
346     .           ( long(ig)*180./pi .ge. 110. ) .and.
347     .           ( long(ig)*180./pi .le. 181. ) )
348     .         .or.
349
350     .         ( ( lati(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
351     .           ( lati(ig)*180./pi .le. 76.  ) .and.
352     .           ( long(ig)*180./pi .ge. 150. ) .and.
353     .           ( long(ig)*180./pi .le. 168. ) )
354     .         .or.
355     .         ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
356     .           ( lati(ig)*180./pi .le. 80.  ) .and.
357     .           ( long(ig)*180./pi .ge. -150.) .and.
358     .           ( long(ig)*180./pi .le. -110.) ) )
359     .         then
360             
361               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
362              !    watercaptag(ig)=.true.
363                  alternate = 1
364               else
365                  alternate = 0
366               endif !end if alternate = 0
367               
368       endif
369
370c----------- Opposite olympia planitia water cap --------
371c--------------------------------------------------------
372
373        if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
374     .         ( lati(ig)*180./pi     .le.  84 ) )
375     .         .and.
376     .       ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
377     .         ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
378!     .     ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
379!     .         ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
380!     .       ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
381!     .         ( long(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
382        !   watercaptag(ig)=.true.
383        endif
384
385
386c -------------------- Central cap ----------------------
387c--------------------------------------------------------
388
389      if (abs(lati(ig)*180./pi).gt.80)
390     .       watercaptag(ig)=.true.
391           
392c--------------------------------------------------------
393c--------------------------------------------------------
394      end do ! of (ngridmx)
395
396
397       ELSE
398     
399         print*, 'In surfini.F, icelocationmode is ', icelocationmode
400         print*, 'It should be 1, 2 or 3.'
401         call abort
402
403       ENDIF ! of if (icelocation)
404       
405       
406        ! print caps locations - useful for plots too
407        print*,'latitude | longitude | ig'
408        do ig=1,ngridmx
409          dryness (ig) = icedryness
410
411          if (watercaptag(ig)) then
412             print*,'ice water cap', lati(ig)*180./pi,
413     .              long(ig)*180./pi, ig
414          endif
415        enddo
416       
417#endif     
418
419       ENDIF ! (caps & water)
420       
421
422c ===============================================================
423c      INITIAL ALBEDO
424c ===============================================================
425
426         write(*,*)"surfini: water frost thickness",
427     s     frost_albedo_threshold
428         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
429         write(*,*)"surfini: water ice TI:", inert_h2o_ice
430
431c        To start with : Initial albedo = observed dataset
432c        -------------------------------------------------
433         DO ig=1,ngrid
434              psolaralb(ig,1)=albedodat(ig)
435              psolaralb(ig,2)=albedodat(ig)
436         END DO
437         PRINT*,'minimum albedo sans water caps',
438     s     albedodat(ISMIN(ngrid,albedodat,1))
439         PRINT*,'maximum albedo sans water caps',
440     s     albedodat(ISMAX(ngrid,albedodat,1))
441
442c        initial albedo if permanent H2O ice is present
443c        ------------------------------------------------
444         IF ((water) .and. (caps)) THEN
445           DO ig=1,ngrid
446            IF (watercaptag(ig)) THEN
447              psolaralb(ig,1) = albedo_h2o_ice
448              psolaralb(ig,2) = albedo_h2o_ice
449            ENDIF
450           END DO
451           PRINT*,'minimum albedo avec water caps',
452     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
453           PRINT*,'maximum albedo avec water caps',
454     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
455
456         ENDIF
457
458c      changing initial albedo if CO2 ice is present
459c      -------------------------------------------
460
461       DO ig=1,ngrid
462         IF (piceco2(ig) .GT. 0.) THEN
463             IF(ig.GT.ngrid/2+1) THEN
464                icap=2
465             ELSE
466                icap=1
467             ENDIF
468             psolaralb(ig,1) = albedice(icap)
469             psolaralb(ig,2) = albedice(icap)   
470         END IF
471       END DO
472
473c      changing initial albedo if water ice frost is present
474c      -------------------------------------------
475       IF (water) THEN
476          do iq=1,nqmx
477c          if there is frost and surface albedo is set to albedo_h2o_ice
478           if(noms(iq).eq."h2o_ice") then
479             do ig=1,ngrid
480             
481              if ((watercaptag(ig).eqv..false.)
482     &     .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
483              print*, ''
484              print*, '!!! PROBLEM in SURFINI !!!!'
485              print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
486     & WATERCAPTAG IS FALSE'
487              print*, ''
488              print*, 'ig,qsurf,threshold' ,
489     &         ig, qsurf(ig,iq), -frost_albedo_threshold
490              print*, ''
491              print*, '1) Check h2o_ice in startfi and ice
492     & distribution in surfini'
493              print*, '2) Use ini_h2osurf option in newstart'
494              print*, ''
495              CALL ABORT
496            endif
497             
498              if ((piceco2(ig) .eq. 0.).and.
499     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
500                     psolaralb(ig,1) = albedo_h2o_ice
501                     psolaralb(ig,2) = albedo_h2o_ice
502c                     PRINT*,'surfini.F frost',
503c     &                  lati(ig)*180./pi, long(ig)*180./pi
504               endif
505              enddo
506           endif
507          end do
508          PRINT*,'minimum albedo avec givre et co2',
509     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
510          PRINT*,'maximum albedo avec givre et co2',
511     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
512       END IF
513         
514      RETURN
515      END
Note: See TracBrowser for help on using the repository browser.